2.1 Control problem

LDR brachytherapy is one of the most effective treatment modalities for early-stage prostate cancer. The treatment result greatly depends on the precise placement of the needle and accurate delivery of radioactive seeds and the dose distribution. However, the anatomy of the prostate, which has a centrally located urethra, demands that the needle should follow an active curvilinear path in order to avoid puncturing the urethra and reducing the dose coverage to the urethra (Podder et al. 2012, Ruiz et al. 2012). As shown in Fig. 1, the entire process required to reach the target point in the world coordinate frame can be divided into two sub-processes, namely the robot external navigation portion and the flexible-needle, internal-target-entering portion. In the external portion, the robot moves the needle tip from the home position point (X₀, Y₀, Z₀) to the point on the body surface (X_i, Y_i, Z_i) with no needle deflection. The 3DOF robot holding the needle is capable of moving in 3D space when outside the patient’s body. In the internal portion, that is, inside the body, the flexible needle is curved by heating the SMA along with the insertion into the body by the robot in order to reach the final target point. The needle tip is moved/bent in the y-direction by actuating the SMA actuator/wire. A planar, curved trajectory with the initial (X_i, Y_i, Z_i) on the patient’s body surface and the final (X_f, Y_f, Z_f) target position are specified to meet the clinical requirement, that is, to perform the desired tasks (delivery of therapeutic agent or collection of biopsy sample, etc.) accurately (Podder et al. 2014). To reach the target in the prostate, 1DOF from the needle bending is moving the needle tip in the y-direction and 1DOF from the robot moving in the x-direction complete the 2D movement, and which can be achieved in a differential velocity mode of the robot and the SMA actuated needle bending. This planar trajectory is the subset of the 3D trajectory which covers each slice of the whole volume of the prostate gland.

Figure 1.

Robot-needle system concept for surgical intervention.

2.2 Kinematics model of the robot-needle system

When the robot-needle system functions in the external portion, that is, moving from (X₀,Y₀,Z₀) to (X_i, Y_i, Z_i), its kinematics are relatively trivial as all of the three actuators of the robot are linear. This kind of 3D robot is termed a Cartesian robot. Its transformation is given as:

(1) \begin{align}{g_{cart}} = \left[ \begin{matrix} 1 &\quad 0 &\quad 0 & \quad {{P_x}}\\[6pt] 0 & \quad 1 & \quad 0 & \quad {{P_y}} \\[6pt] 0 & \quad 0 & \quad 1 & \quad {{P_z}} \\[6pt] 0 & \quad 0 & \quad 0 & \quad 1 \end{matrix} \right]\end{align}

where P _x, P _y, and P _z are translations in the three directions. As the robot-needle system begins to enter the human body, the active flexible needle has two DOFs, that is, penetration into the body and bending of the tip, both of which are actuated by a robot joint motor angle in the X-axis and the input voltage to the SMA, respectively.

The robot continues to push the needle into the tissue as long as the voltage supplied to the SMA causes the needle tip bend and thus reaches the target point in a desired curvilinear trajectory. The different radii of the paths caused by various input voltages are attained using different bending of the tip. We assume that the needle tip is bending within the x-y plane. The value of ${\alpha _0}$ indicating the initial value of the needle tip orientation off the x-y plane should be calculated.

(2) \begin{align}{ \propto _0} = \left\{ {\begin{matrix} {\arctan \left( {{\raise0.7ex\hbox{${{y_n}}$} \!\mathord{\left/ {\vphantom {{{y_n}} {{z_n}}}}\right.}\!\lower0.7ex\hbox{${{z_n}}$}}} \right),{z_n} \ge 0} \\[6pt] {{\rm{\pi }} + \arctan \left( {{\raise0.7ex\hbox{${{y_n}}$} \!\mathord{\left/ {\vphantom {{{y_n}} {{z_n}}}}\right.}\!\lower0.7ex\hbox{${{z_n}}$}}} \right),{z_n} \lt 0} \\ \end{matrix} } \right.\end{align}

Figure 2.

Needle curving coordination.

The position of the needle tip could be viewed as the rotation of the virtual axis according to the radius and angle, as shown in Fig. 2. If we divide the needle-tip path into a series of segments,${C_1}, \ldots ..,{C_n}$, we have parameterized ${U_i}\left( {{r_i},{l_i},{\phi _i}} \right)$ in each segment of ${C_i}$, where r _i is the radius of the i ^th arc (corresponding to input voltage), l _i is the length of the i ^th segment (corresponding to the robot placement in the x-axis), and Φ is the degree of the arc. The transformation matrix for each g _i segment can be formulated as follows:

(3) \begin{align}{g_i} = \left[ {\begin{matrix} {\cos {\phi _i}} & \quad 0 & \quad {\sin {\phi _i}} & \quad {{x_i}} \\[6pt] 0 &\quad 1 &\quad 0 &\quad 0 \\[6pt] { - \sin {\phi _i}} &\quad 0 &\quad {\cos {\phi _i}} &\quad 0 \\[6pt] 0 &\quad 0 &\quad 0 &\quad 1 \end{matrix} } \right]\end{align}

where ${{\phi _i}}$ is the angular position between two adjacent segments. The entire transformation from the input joint of the robot to the needle-tip position is given by:

(4) \begin{align}g_j^t = {g_{cart}}\mathop \prod \limits_{i = 1}^{i = n} {g_i}\end{align}

The inverse kinematics of the external robot portion is the inverse of the direct kinematics as it is linearly decoupled. After the derivation of the ${\alpha _0}$, the robot-needle inverse kinematics could be formalized by the geometric relationship shown in Fig. 2, where different parameters are defined below.

(5) \begin{equation}r_i=y_n/2+(x_n^2)/2\end{equation}

(6) \begin{equation}\varnothing_{i}=\pi-2arctan(x_{n}/2)\end{equation}

(7) \begin{align}{l_i} = {r_i} \times {\emptyset _i}\end{align}

2.3 Control method

Figure 3 demonstrates the control strategy for the two-phase needle insertion. As there is no time overlap between the two phases, the x-axis robot could be used for both the internal and external movement. We divided the control plan into two subsections. In the first one shown in the red dash line area, the 3D Cartesian robot would move the needle tip from the home position (X₀, Y₀, Z₀) to the initial insertion position (X_i, Y_i, Z_i) on the patient’s skin. Three PID controllers would separately serve the three control loops. In the second subsection, we introduce the coordinated control strategy indicated in the green dash area and involving a two-loop control strategy in which the needle control loop regulates the deflection of the SMA-actuated flexible needle.

Figure 3.

Robot-needle system control diagram.

The robot insertion loop regulates the robot guidance by pushing the needle along the x-axis, while the SMA actuator is bending the needle to move the needle tip in the y-direction. These two combined x-direction and y-direction motions of the needle tip will make the needle follow the desired curved trajectory inside a patient’s body. The curvilinear relationship between the robot position (${R_x}$) and the needle position (${N_y}$) is defined by the relationship demonstrated below:

(8) \begin{align}{R_x} = \propto {N_y}\end{align}

The measured position executed by the robot would provide the input to the needle loop in order to generate the curving. The EM sensor would provide measurement of the tip position in each axis. Therefore, the computed positional errors, that is, the difference between the desired and measured positions of the needle tip, are provided to the corresponding controller. Hence, the reference and output relationship of the coordinated control system in the green dash-line box, that is, the x- and y-direction movements, would follow the desired trajectory.

2.3.1 Robot-steering part: 3D decoupled control in Cartesian coordination

The robot part of the robot-needle system is Cartesian whose three principal axes of control are linear and are at right angles to each other. Three motors drive the three sliding joints of the robot to move the needle-tip position in a straight line in the wrist up-down (y), in-out (z), and left-right (x).

The robot is controlled using the digital PID controller whose equations are given in ((9)).

(9) \begin{align}{u_{PID}}(k) = {u_{PID}}( {k - 1}) + {K_P}\left[ {e(k) - e( {k - 1})} \right] + {K_I}{T_S}e(k)\,{ + _D}\left[ {{e_f}(k) - 2{e_f}( {k - 1}) + {e_f}\left( {k - 2} \right)} \right]\end{align}

where ${T_S}$ is the sampling period and $e(k) = {x_d}(k) - x(k)$ is the error signal at the ${k^{th}}$instant. ${e_f}(k)$ is the first-order, low-pass-filtered error signal given by ${e_f}(k) = \frac{{{\tau _f}}}{{{\tau _f} + {T_S}}}{e_f}( {k - 1}) + \frac{{{T_S}}}{{{\tau _f} + {T_S}}}e(k)$, ${\tau _f}$ is the time constant of the filter, and $u(k)$ is the controller output at the ${k^{th}}$instant. ${K_P}$, ${K_I}$, ${K_D}$ are the PID gains, respectively.

2.3.2 Flexible-needle part: 1D robot actuation with 1D needle curving

1. A. System description

The dynamic model of the SMA-actuated, flexible needle consists of a spring-mass-damper system (Hill et al., 2009). Its equation is:

(10) \begin{align} \epsilon = \frac{f(t)-k_{s}L_{\epsilon}-cL\epsilon-A\sigma}{mL}\end{align}

where $\epsilon=x/L$ is the strain in the needle actuated by ${\varepsilon _s}$ the strain of the SMA wire. $A$ and $\sigma $ are the area and stress of the SMA wire, respectively. ${k_s}$ is the spring constant, $c$ is the damping constant, $m$is the mass, and $f(t)$is the externally applied force.

From the control system perspective, for the flexible needle, we have the voltage (current equivalent) as input to drive the needle-tip position, which is the output. The mechanism of the SMA actuation is governed by the heat generated by the input voltage to the SMA wire. The applied voltage and change in the SMA wire temperature are represented by the equation below (Songye and Yunfeng, 2007).

(11) \begin{align}\dot T = \frac{{{k_h}}}{{A{p_0}{C_p}{L_{SMA}}}}\left( {T - {T_0}} \right) + \frac{{{V^2}}}{{A{p_0}{C_p}R}}\end{align}

where ${k_h}$, ${p_o}$, ${C_p}$ are the heat convection coefficient, the density, and the specific heat of the material, respectively. ${T_0}$ is the ambient temperature, $V$ is the applied voltage, and $R$ is the electrical resistance of the SMA wire.

Due to the change in the SMA wire temperature, the induced stress of the wire also changes and which is approximated by the thermo-mechanics relationship presented below (Prahlad and Chopra, 2003).

(12) \begin{align}\dot \sigma = E\left( \delta \right){\dot \epsilon _s} - { \epsilon _1}E\left( \delta \right){\dot \epsilon _s} + \left( {{E_m} + {E_a}} \right)\left( {{ \epsilon _s} - { \epsilon _1}\delta } \right){\dot \epsilon _s} + \theta \dot T\end{align}

where $E$ is the Young’s modulus of the SMA, $\delta $ is the martensite fraction, $\theta $ is the thermal coefficient, and ${\varepsilon _1}$ is the maximum residual strain. ${E_a}$ and ${E_m}$ are the Young’s modulus for austenite and martensite, respectively.

The stress would cause the strain of the SMA wire. Finally, the needle-tip position changes along the strain of the wire due to the internal physical effects. Combining ((11)) and (12) in addition to differentiating (10) derives a more general form as follows:

(13) \begin{align} \dddot{\epsilon}= d\left( {\sigma ,T, \epsilon ,\dot \epsilon ,\ddot \epsilon } \right) + \dot f(t) + b\left( {\sigma ,T, \epsilon ,\dot \epsilon ,\ddot \epsilon } \right)u\end{align}

where $u$ is the control input. We could transform the system model (13) into the state space representation form in terms of states of the flexible needle, as follows:

(14) \begin{gather}{\dot \chi _1}(t) = {\chi _2}(t)\nonumber\\[3pt] {\dot \chi _2}(t) = {\chi _3}(t)\nonumber\\[3pt] {\dot \chi _3}(t) = d(\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t) + \Delta d(\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t) + b(\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t)u + f(t)\nonumber\\[3pt] y = {\chi _1}(t)\end{gather}

where ${\chi _1}(t)\ {\chi _2}(t),$and ${\chi _3}(t)$ are the system states representing the position, velocity, and acceleration of the SMA actuator, respectively; $u$ and $y$ are the control input and system output, respectively; $b\left( {\sigma ,T,{\chi _1},{\chi _2},{\chi _3}} \right)$ is the input gain; $b\left( {\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t} \right)$ and $\Delta d\left( {\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t} \right)$ are the system nominal parameter and uncertainty, respectively. $f(t)$ represents the external disturbance or loading. The assumption is that $\Delta d\left( {\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t} \right)$ and $f(t)$ are bounded by two parameters, α and ${d_0},$ such that the conditions, $\left| {\Delta d(\sigma ,T,{\chi _1},{\chi _2},{\chi _3},t)} \right| \le {d_0}$ and $\left| {f(t)} \right| \le \propto $, are satisfied. It is also assumed that the input gain, $b\left( {\sigma ,T,{\chi _1},{\chi _2},{\chi _3}} \right)$ is certain to be within the range of ${b_{min}} \le b \le {b_{max}}$.

1. B. Controller design

   1. a. Proportional, integral, and derivative with the nonlinearity (PID-P³) controller

We deployed the PID-P³ controller to deal with the nonlinearity of the flexible-needle position control. The discretized PID-P³ control equation used in Maria Joseph et al. (2014) is given by

(15) \begin{align}{u_{pid - p3}}(k) &= {u_{pid - p3}}( {k - 1}) + {K_P}\left[ {e(k) - e( {k - 1})} \right] + {{{K_P}h} \over {{T_I}}}e(k) + {K_P}{T_D}[{e_f}(k)\,\nonumber\\ &\quad- 2{e_f}( {k - 1}) + {e_f}\left( {k - 2} \right)] + {K_{p3}}{\left[ {e(k) - e( {k - 1})} \right]^3}\end{align}

where h is the time step; $e(k) = {x_d}(k) - x(k)$ is the error at the $k$th instant; ${e_f}(k)$ is the low-pass filtered error given by ${e_f}(k) = \frac{\tau }{{{\tau _f} + {T_S}}}{e_f}( {k - 1}) + \frac{{{T_S}}}{{{\tau _f} + {T_S}}}e(k)$; $u(k)$ is the controller output at $k$th instant; and ${K_P}$, $\frac{{{K_P}h}}{{{T_I}}}$, ${K_P}{T_D}$ are the PID gains, respectively.

1. b. Static PID Sliding Mode Controller (SPID-SMC)

Sliding Mode Control

In this study, we designed a position controller for the flexible needle based on the Sliding Mode theory. From Slotine and Li (1991), a sliding surface for a third-order system can be defined as:

(16) \begin{align}s = \ddot \chi + 2\mu \dot{\tilde{z}} + {\mu ^2}\tilde \chi \end{align}

where $\tilde z = {z_d} - z$, is the desired position of the flexible needle tip. $\mu $ is a strictly positive constant. The equivalent control is defined by stating that the system will not leave the sliding surface once it reaches the surface.

(17) \begin{align}\dot s = 0\end{align}

By differentiating the equation and applying the sliding condition, we can derive the equivalent control law by approximating the dynamics.

(18) \begin{align}{u_{eq}} = {\hat b^{ - 1}}\left[{\dddot{\chi}_d - d - \nabla d - f + 2\mu \ddot \chi + {\mu ^2}\dot{\tilde{\chi}}} \right]\end{align}

The equivalent control law allows the system states to be on the sliding surface. On the other hand, we still need the discontinuous control, ${u_N}(t)$, which handles the system uncertainties, especially for the dominating nonlinear hysteresis, and guarantees the stability of the system by switching the system to the sliding mode surface as well.

(19) \begin{align}{u_N} = - {K_{SMC}}\ {{\ast}}\ sat(s)\end{align}

where

(20) \begin{align}sat(s) = \left\{ {\begin{matrix} 1 &\quad {s \gt 1} \\[6pt] { - 1} &\quad {s \lt - 1} \\[6pt] s &{others} \end{matrix} } \right.\end{align}

In general, the total control input is given as follows:

(21) \begin{align}u = {u_N} + {u_{eq}}\end{align}

The dynamics model of the system is too complicated to describe accurately enough via a set of the nonlinear differential equations. An efficient PID equivalent controller is shown below.

Static PID approximation

The approximate unknown and uncertainty parameters in the equivalent control are achieved using a PID controller; here we call it static PID (SPID) compared to the adaptive one. In (22), ${u_{eq}}$ also includes the derivative and proportional term regarding the error of the needle position. Here we present the integral element in order to increase the control performance as the sliding surface converges to zero. Therefore, we give the following equation:

(22) \begin{align}\eqalign{ {u_{eq}} \approx {u_{APID}} = {K_P}{\rm{e}} + {K_I}\mathop \int e\,dt + {K_D}{{\dot {\textrm{e}}}} \cr {{\rm{u}}_{SPID - SMC}} = {u_{SPID}} + {u_N} \cr} \end{align}

1. c. Robust Adaptive PID Sliding Mode Controller

The adaptive PID has an added advantage over the static one. It can well undertake the equivalent control by generating the control signal even in the case of hysteresis and disturbance; it is ensured that the sliding surface goes to zero, and thus the chattering problem prevalent in the sliding mode control is suppressed. If we apply the adaptive laws for the PID control gains, ${K_P}$, ${K_I}$, ${K_D}$ may be determined using the gradient descent method:

(23) \begin{align}{K_P}( {k + 1}) = {K_P}(k) + {\gamma _1}se\nonumber\\[3pt] {K_I}( {k + 1}) = {K_I}(k) + {\gamma _2}s\int edt \nonumber\\[3pt] {K_D}( {k + 1}) = {K_D}(k) + {\gamma _3}s\dot e \end{align}

in which ${\gamma _1},\,{\gamma _2},\,{\gamma _3}$ are the learning rates. As ${u_{APID}}$ becomes equal to ${u_{eq}}$, $s \to 0$. It is well known that the gradient parameter adaption is not sufficiently robust when the system is subject to bounded disturbances. We make some modifications of the adaptive controller in (24) by adding the term $ - {\sigma _s}K(k)$ to improve robustness, where ${\sigma _s}$ is a positive constant parameter and the gainis the vector form for ${K_P}$, ${K_I}$, ${K_D}$ (see (24)). The is defined as the switching condition introduced in order to achieve robust asymptotic tracking performance (Ge et al., 1998), where ${\left\| {\rm{K}} \right\|_{\max }}$ is the maximum value of $\left\| {\rm{K}} \right\|$. In (25), we could present the control law of RAPID-SMC.

(24) \begin{gather}{K_P}(k) = \left( {1 - {\sigma _s}} \right){K_p}( {k - 1}) + {\gamma _1}se\nonumber\\[3pt] {K_I}(k) = \left( {1 - {\sigma _s}} \right){K_I}( {k - 1}) + {\gamma _2}s\smallint edt \nonumber\\[3pt] {K_D}(k) = \left( {1 - {\sigma _s}} \right){K_D}( {k - 1}) + {\gamma _3}s\dot e\nonumber\\[3pt]{\sigma _s} = \left\{ {\begin{matrix} {0,||K|| < ||K|{|_{max}}} \cr {{\sigma _0}\left( {{\raise0.7ex\hbox{${||K||}$} \!\mathord{\left/ {\vphantom {{||K||} {||K|{|_{max}} - 1}}}\right.}\!\lower0.7ex\hbox{${||K|{|_{max}} - 1}$}}} \right),\,||K|{|_{max}} \le ||K|| \le 2||K|{|_{max}}} \cr {{\sigma _0},||K|| \ge 2||K|{|_{max}}} \end{matrix} } \right. \end{gather}

(25) \begin{align} {u_{eq}} \approx {u_{APID}} &= {K_P}(k){\rm{e}} + {K_I}(k)\mathop \smallint \nolimits^ e{\rm{}}dt + {K_D}(k){\dot{\textrm{e}}}\nonumber\\ &\quad {{\rm{u}}_{APID - SMC}} = {u_{APID}} + {u_N} \end{align}

3.1 Experiment setup

The experimental setup for the 3D robot with an SMA-actuated, flexible needle is shown in Fig. 4. A 5DOF EM tracking sensor (Aurora, Northern Digital Inc., Waterloo, Ontario, Canada) measures the actual needle-tip position at a sampling period of 1 ms (Fig. 5) and depicts the zoomed-in view of the SMA wire actuated by the active flexible needle prototype used in our experiment.

Figure 4.

Experimental setup.

Figure 5.

Zoomed View of the SMA-actuated, active flexible needle.

The overall experimental system parameters are listed in Table I.

Table I.

List of experimental parameters.

The EM sensor, insulated coil is attached to the needle tip. The zero point of the global reference frame is at the top center of the EM field generator. The interface between the active flexible needle’s power amplifier circuit and the control computer is the QUANSER card. We implement the proposed algorithms in the LabVIEW (National Instruments Corporation, Dublin, Ireland) programming environment on the Dell XPS control computer, as described in Table I. We simulate two setups for robot motion only from the home position to the initial insertion position and then to reach the target position through two, different planar motion trajectories for the robot needle system to follow during the insertion phase. In clinical practice, it is desirable to have a well-defined home position of the robot along with the surgical needle. This home position must be away from the patient and at a safe position in the clinical workspace. The two, path settings with each of the two, distinguished motion phases are listed at the bottom of Table I. The three axes of the robot work simutanously delivering the needle-tip position in the first phase, that is, from the home position to the position where the tip of the needle just touches the patient’s skin at the desired location. The second phase of each insertion into the soft-tissue phantom lasts for 5 s.

3.2 Comparative experimental results and discussion

We have implemented the control methods of the robot needle system in two, different path settings for the purpose of its evaluation. In the first phase of each path setting, the robot drives the needle tip from the home position to the initial position of insertion. During this phase, the PID controller is deployed. In the second phase of each path setting, the robot and the flexible needle move in a coordinated way so that the needle tip moves from the initial insertion point to the final target point by following the given trajectory. Although all of the control algorithms in this phase are non-model based, we still investigate their potential properties as we review at the dynamics of the flexible needle. We characterize the four, control methods according to their unique features: (1) without considering the system nonlinearity and without adaptive compensation, namely PID; (2) considering the system nonlinearity but without uncertainty compensation, namely PID-P³; (3) considering both the system nonlinearity and uncertainty compensation but with no adaptive mechanism, namely SPID-SMC; and (4) considering both nonlinearity and uncertainty compensation as well as applying a robust adaptive mechanism, namely RAPID-SMC.

We have attuned the parameters in all of the control algorithms so as to minimize the RMS errors, using the two kinds of needle-tip, curvilinear reference trajectories. We summarize the parameters of the algorithms in Table II. We also give some guidelines regarding tuning the principle parameters. For the robot actuation phase, the PID controller parameters are first determined. We choose the ${K_{pr}}$ and ${K_{ir}}$ to be the primary parameters and keep the ${K_{dr}}$ very small (sometimes at zero) in order to minimize the effects induced by high-frequency noises for both the robot loop and the needle loop. For the control parameters in the needle control loop, PID, and PID-P³, due to the nonlinearity of the needle, we determine the ${K_{pn}}$, ${K_{in}}$, and ${K_{p3}}$around the point where the needle tip begins to rise up curving from a relatively flat plateau and is then corrected using the trial and error method. For the SPID-SMC controller, we give PID most of the control authority and less to the SMC portion first and correct the parameters using the trial and error method. Whereas for the RAPID-SMC, we select ${K_{pn}}$, ${K_{in}}$ and the learning rate of each of the gain parameters of P, I, and D less than 0.5, such that it would keep the system stable, and then attune the main parameters of the SMC controller in order to fill the little gap left by the RAPID controller, and thus achieving the desired reference trajectory.

Table II.

Parameters of the controller in the robot needle system.

The robot needle system targeting process includes two phases. The first phase is the robot driving the needle tip in air. The second phase is the flexible-needle insertion positioning in a tissue phantom. We summarize the two phases of the RMS tracking control results and the algorithm variance analysis in Tables III and IV, respectively. The RMSE position of the tracking control and the RMS of the control action, which is voltage applied to the SMA or motor, with a work length of $n$ are given by (26) and (27), respectively.

(26) \begin{align}{\rm{RMSE}} = \sqrt {\frac{1}{{\rm{n}}}[\left( {\sum {{\rm{x}}_{{\rm{err}}1}}^2 + \sum {{\rm{y}}_{{\rm{err}}1}}^2 + \sum {{\rm{z}}_{{\rm{err}}1}}^2} \right)} \end{align}

(27) \begin{align}{\rm{RM}}{{\rm{S}}_{{\rm{robot}}}} = \sqrt {{1 \over {\rm{n}}}\sum control\ action_r^2} ,\,\,{\rm{RM}}{{\rm{S}}_{{\rm{needle}}}} = \sqrt {{1 \over {\rm{n}}}\sum control\ action_n^2}\\[-20pt] \nonumber \end{align}

Table III.

Robot-delivery PID tracking control results in air.

Table IV.

Needle insertion tracking the control results in the phantom.

The “control action” in (27) refers to the input voltage to the robot and/or the voltage supplied to the SMA actuator. The entire working scenario of the robot needle system under the first path setting using the RAPID-SMC algorithm is shown in Fig. 6. The different phase-motion tracking results are shown in Table V where we separately present the comparative curves in air (Fig. 7) and in the phantom (Fig. 8).

Figure 6.

Robot needle system motion under the first path setting.

Figure 7.

(a) Robot delivery motion in the x-axis; (b) Robot-needle positioning in the y-axis; (c) Robot delivery motion in the z-axis; (d) Robot delivery 3D motion profile.

Figure 8.

(a) Robot-needle positioning in the x-axis; (b) Robot-needle positioning in the y-axis; (c) Tracking control results of RAPID-SMC with PID; (d) Tracking control results of RAPID-SMC with PID-P$^{3}$.

It should be noted that motion control results of the robot in the air delivery shown in Table III and Fig. 7 under both path settings do not have large differences in terms of the trajectory tracking errors. The performance of the second phase shown in Table IV under each path setting differentiates each other as the needle insertion errors in the phantom are much larger than those in the robot delivery phase in air. However, the needle and robot motion in air are clinically very important in order to move the needle from a safe home position to the desired contact position on the patient. We will now increasingly focus on the discussion results of the in-phantom needle insertion experiment and which are more challenging.

According to the results shown in Table IV, in general, the RMSE of the second path setting is larger than those of the first one. The greater curvature of the second path than of the first path was more challenging for the flexible needle to track due to its greater nonlinearity. Comparing all of the four control algorithms, the RAPID-SMC achieved the smallest tracking RMS errors and had the best control action efficiency. The RAPID-SMC has a quicker and more accurate response to the time-varying change corresponding to the structural nonlinearity, due to its hysteresis compared to other algorithms. Moreover, RAPID-SMC is also more capable of handling a situation where the prior knowledge of the needle dynamical model is insufficient. As far as we know, the main difficulties in the control process are due to the nonlinear dynamics of the flexible needle as well as the uncertainties induced by both the deflection of the needle and deformation of interaction with tissue. The reason the RAPID-SMC could outperform all of the other listed algorithms is due to its own two, inherited features, that is, nonlinearity-dealing property and uncertainty-dealing property. We will now discuss the impact that these features have on the experimental results.

1. (1) Nonlinearity property

   • Structural nonlinearity

     The flexible needle has inherent structural nonlinearity which can be represented by a series of nonlinear differential equations.

   • Hysteresis variation nonlinearity

Other than PID, all of the other three control methods can compensate for the nonlinearity changes in different ways. The PID-P³ could balance itself via the cubic nonlinear term. SPID-SMC and RAPID-SMC use the slide-mode techniques maintaining the error and its derivative on the slide surface in order to realize the nonlinear, dynamic structure switching for hysteresis compensation. What causes the difference between the control results is the setup of the matched mechanism in each algorithm. As expected, the PID has a linear tool that could not successfully function with the nonlinearity. The PID-P³ has a slightly nonlinear handling capability, while the SPID-SMC and RAPID-SMC have more powerful solutions in order to take into account fast time-varying changes in the needle-tip position. Figure 8(c) demonstrates the difference in the position-tracking error between PID and RAPID-SMC using the first trajectory. The SPID-SMC and RAPID-SMC overcome the shortcomings of the PID and PID-P³ in both curvilinear trajectory tracking methods. The discontinuous control term in SMC is mainly used to manage the more demanding curving situation, especially for achieving the highly curved, clinical configurations. Adding the discontinuous term provides the mechanism to conquer the hysteresis behavior such that SMC-related methods can achieve accurate control when the reference trajectory changes rapidly. From Table IV and Fig. 8, we observe that the control errors using the nonlinear control method (SPID-SMC and RAPID-SMC) are approximately 27% less than those of the linear algorithm (PID).

1. (2) Uncertainty dealing property

   • Model uncertainty

When we face the inadequate prior model knowledge of the SMA-actuated, flexible needle and we have to manage the various time-dependent changes, the uncertainty-dealing property is required. Performance of the RAPID-SMC is superior compared to the other algorithms. RAPID-SMC puts forth an adaptation control strategy in order to compensate for the deficiencies of the non-model by adaptively changing the PID parameters. In contrast, all of the other algorithms lack the mechanism to take into account the unexpected changes and uncertainties. The tracking control results of RAPID-SMC in Table IV from different arc trajectories and in Fig. 8(d) show the robustness of the RAPID-SMC compared with PID-P³ in more detail together with the smaller time scale. In Fig. 8(c), we observe that RAPID-SMC made the needle preserve the same level deviation when it encountered the sharper path/trajectory, whereas the PID-P³ became worse and could not avoid a certain degree of vibration in the needle tip.

• Tissue interaction uncertainty

As the needle moves further into the tissue phantom, the interactive disturbance between the needle and the tissue environment is more severe. The errors increase with the depth of the needle in the tissue phantom (Fig. 8(d)). RAPID-SMC is better for handling the disturbance loading as the adaptive mechanism could continuously work against the disturbance and thus minimize the errors within the stable boundary. Other controllers are not capable of dealing with the issue of interaction disturbance.

We also found that the high accuracy of the robot loop control could improve the accuracy of the needle loop control. The smooth robot position output could induce a stable trajectory tracking; however, the noisy robot position output introduced more errors in the insertion process. In certain even worse cases, the noisy robot position could result in vibration of the needle tip positioning and unstable control behavior. For example, when the robot drove the needle tip to enter the phantom from air to phantom, a resisting force in the opposite direction against the needle insertion would pose a difficulty for the PID controller, and thus introduces larger errors. Therefore, the highly stable and high precision control of the robot loop is important in the coordinated control.