# Prediction Control for Brachytherapy Robotic System

Authors: Ivan Buzurovic; Tarun K. Podder; Yan Yu
Journal: Journal of Robotics
 (2010)

---

## Abstract

In contemporary brachytherapy procedure, needle placement at desired location is challenging due to a variety of reasons. We have designed and fabricated an image-guided robot-assisted brachytherapy system to improve the needle placement and seed delivery. In this article we have used two different predictive control strategies in order to investigate the needle insertion efficacy and system dynamics during prostate brachytherapy. First, we used neural network predictive control (NNPC) to predict an insertion force. The NNPC uses the linearized state-space model of the robotic system to predict future system performances. Second, we used feedforward model predictive control (MPC) which allows the controller to compensate the influence of a measured disturbance's impact immediately rather than waiting until the effect appears in the system. Feedback control problem for the contact force regulation is considered.  The simulation results and experiments for both cases are presented and compared.


---

## Body

## 1. Introduction

LOW-DOSE rate (LDR) prostate brachytherapy is a method of delivering radiotherapy by implanting radioactive sources into and around the prostate gland [1]. This procedure is performed under the guidance of transrectal ultrasound (TRUS) images. In traditional brachytherapy procedures, the needles are inserted transperineally under the guidance of transrectal ultrasound images. Both the needle and the ultrasound are operated manually. The seeds are deposited using a manual applicator.  In order to increase accuracy of the needle placement and seed delivery, we have developed a fully automated robotic system for prostate brachytherapy [2]. Using the robotic approach we are able to record the needle insertion forces and motion trajectories measured during the actual brachytherapy needle insertion while implanting radioactive seeds in the prostate gland [3, 4]. These insertion forces are significantly responsible for needle deviation from the desired trajectories and target movements. Also, needle insertion causes soft tissues to deform and it is often difficult to obtain precise imaging data during insertion. Since the target organ is deformable, sometimes it is necessary for the next needle trajectory to be updated according to the magnitude of the deformation. During brachytherapy needle insertion forces deform the prostate tissue and displace the targeted seed positions [5].  The prostate deformation and possible calcification within the gland can cause additional needle deflection. Several authors discussed these problems and suggested methods and procedures to improve the needle insertion efficacy.In order to reduce tissue deformation, the authors have studied the effect of different trajectories for a 2-DOF robot performing needle insertion in soft tissue [6]. They have compared tissue deformation and infinitesimal force per tissue displacement for different trajectories. In the articles [7, 8] the needle insertion motion planning system is presented. The system is based on an interactive simulation of needle insertion in deformable tissues and numerical optimization to reduce placement error. The authors described a 2D physically based, dynamic simulation of needle insertion that uses a finite-element model of deformable soft tissues and models needle cutting and frictional forces along the needle shaft. The investigation of flexible needle insertion was presented in [9]. A method for an interactive virtual needle insertion simulation has been suggested in [10]. In [11, 12] a force model for needle insertion and a model that predicted the soft tissue deformation caused by needle insertion were presented. A novel method for the evaluation and correction of brachytherapy needle deflection was introduced in [13]. The authors discussed about the influence of needle rotation to the insertion force and needle deformation. It has been reported that robot used in conjunction with an optimized needle insertion technique benefits patients with increased accuracy, leading to a more successful outcome and reduced complications [14].  To improve needle insertion outcome, the needle insertion point, needle heading, and needle depth are optimized by minimizing the distance between a rigid needle and a number of targets in the tissue as it is shown in [15]. In [16] the authors suggested fast needle insertion to minimize tissue deformation and damage.Generally, to improve needle placement and seed deposition in brachytherapy procedure several methods have been presented in the literature, such as parameter optimization, different needle rotation techniques, robotic insertion, force modeling, and needle steering techniques.  In the previous studies, we concluded that the proper selection of the translational and rotational velocities may reduce tissue deformation and target movements by reducing insertion force [17–19]. Therefore, it can be concluded that the insertion force has a dominant influence on the needle insertion, seed deposition precision, and dosimetry distribution in brachytherapy procedure.In our initial work [20] we have investigated tracking problem for the same robotic system using neural network approach. The force prediction model was introduced as well. In this article we have introduced novel methods to control the brachytherapy robot with the lowest tissue deformation and highest seed precision delivery as the ultimate goals. In order to achieve the desired dynamic behavior of the robotic system, both control strategies that we implemented in the robotic system have a force prediction module. For the first one we have used an artificial neural network (ANN) and neural network predictive control (NNPC). The second one was a feedforward model predictive control (MPC). The purpose of both approaches is to control the reactive force which is responsible for tissue displacement. When the reactive force is minimized, tissue displacement decreases. The force prediction control was also used to minimize the effect of system time-delay. The mathematical model of the system should include the contact force. Consequently, the purpose of both NNPC for force prediction and feedforward MPC is to optimize insertion force which should result in less tissue displacement and deformation, on one side, and higher seed deposition precision, on the other side. The experiments and results for both control approaches are discussed in the following parts.In the first approach, we used an NNPC to predict insertion force in order to achieve real-time adaptive needle control. We have tested whether the system is able to predict the insertion force. The implemented optimization algorithm then computes the control signals that optimize the future system performances. It is shown in the literature that the force prediction technique for needle insertion can improve insertion outcome. The measurement and prediction of insertion force and needle fracture force using experimental approach were investigated in [21, 22]. Prediction of how the skin deforms upon insertion by microneedles is described in [23]. Applications of ANN in medical robotics are used in the following studies. Nonetheless, no application of the neural networks in the needle insertion optimization has been reported so far.  ANN is applied in early detection of prostate cancer [24, 25]. In miniature robotic surgical systems ANN is used in conjunction with real-time visual feedback to “learn” the inverse system dynamics and control the manipulator endpoint trajectory [26]. The real-time control of a robot arm using recorded neurons in the motor cortex together with mathematical transformations including neural networks is presented in [27]. The overview of the ANN applications in many disciplines, especially in medicine can be found in [28]. In the article [29] the authors solved the problem of achieving high-accuracy positioning of a medical robot using neural network for the patient positioning system. Another robotic system for cardiac surgery which uses a recurrent neural network with adaptive internal states is described in [30]. In [31] it is shown that the neural network demonstrated robust adaptability to all of the observed breathing patterns while the linear filter failed in a significant percentage of the cases.In the second approach (the MPC) the control strategy was to compare predicted robot states to a set of objectives and to adjust the actuators to achieve the objectives respecting the robot’s constraints. Such constraints included the actuators’ physical limits, needle deflection, tissue deformation and displacement, and the needle contact forces. Therefore, the robotic system is modeled to take the contact with an environment into account. For a described case, robot dynamics is usually called constrained system dynamics. Mathematical models of robotic systems give rise to a mathematical system of differential equations and algebraic equations that can be viewed as a singular system of differential equations [32]. While the manipulator end-effector is in contact with the environment, system dynamics changes significantly.  Notwithstanding, it is essential to maintain an appropriate control over the working regime not only force and torques arising due to contact, but also the position and orientation of the end-effector [33]. In this approach, the controller predicts how much each output will deviate from its setpoint within the prediction horizon. It multiplies each deviation by the output’s weight and computes the weighted sum of squared deviations. The critical step for MPC is choosing the weights. Because of the contact force influence the controller cannot be exclusively focused on setpoint tracking. In that case, the large manipulated-variable adjustment cannot be implemented. Contact forces are capable to accelerate equipment wear or can lead to control system instability. The described problem can be solved by decomposition of the system to its slow and fast subsystems.


## 2. System Description

We developed a 16- (DOF-) robot-assisted brachytherapy system, [2]. This robotic system is divided into three subsystems: the cart, the supporting platform, and the surgery module, with 3, 6, and 7 DOF, respectively, Figure 1. The supporting platform connects the surgery module to the cart. The surgery module consists of a 2-DOF ultrasound probe driver and a 5-DOF needling module. The ultrasound (US) module can be translated and rotated independently by two DC servomotors fitted with high-resolution optical encoders and gearboxes. In this current study, we have investigated the 5-DOF-needling module which consists of a gantry and a needle driver. The gantry connects the needle-driving module to the positioning platform. The gantry has two translation motions and one rotational motion (pitching). The needle driver subsystem consists of a hollow needle (cannula) and a solid needle (stylet) driven separately by two DC servomotors. The cannula rotates continuously or partially using another tiny DC motor. The basic task of these parts is to deliver the exact prescribed dose of radioactive seeds with high precision level into the human prostate. The seeds are delivered through the cannula. During the operation the stylet pushes the seeds through the cannula according to the control algorithm and prescribed surgery plan.Figure 1
Robot-assisted brachytherapy system (Euclidian). Three basic subsystems are shown:  7DOF robot subsystem (surgery module), supporting platform and cart.Also, the system is designed to take ultrasound images during the operation to update the real-time radiation dose distribution, seed position, and number of needles to be inserted into the prostate, depending on the surgery plan. Dedicated software for 3D imaging and control is developed to support the surgery procedure.


## 3. Mathematical Framework

### 3.1. System Dynamics for NNCP

Prior to implementing one of the suggested controllers to the system, it is essential to investigate the system behavior. Therefore, in this section, we have presented a mathematical model of the robotic system together with the mathematical model of the actuators. The derivation of dynamic equations is performed using the energy-based Lagrangian dynamic formulation. The mechanical part of the robotic subsystem is represented as a five-link open kinematic chain. All links are rigid bodies. We assumed that all materials for the system are homogenous. The generalized coordinates (q1, q2, q3, q4, q5) are adopted to describe each joint movement; see Figure 2. The general form of dynamic equation can be written as follows:
(1)∑α=1naαγ(q)q̈α+∑α=1n∑β=1nΓαβ,γ(q)q̇αq̇β=Qγg+Qγu,
where γ=1,…,5, aαγ(q) stands for the metrics tensor coefficients, Γαβ,γ(q) denotes the corresponding coefficients, and Qg, Qu are the generalized forces components. The former is a component due to gravitation and the latter is a component that corresponds to the actuator forces and torques. Each robot joint is motorized by a direct current motor, depending on the torque. Mathematical relationship between the output torque and the input motor torque for each motor is(2)NVNmJMθ̈+Fθ̇+M=CnNmIR.
Equation of the electrical equilibrium in the rotor circuit is: (3)LRİR+RRIR+CENVθ̇=U,
where θis the output rotation angle, Mis the output torque of the motor, IRis the rotor current, LRis the motor’s inductance, and U is the motor’s voltage. Other values are constant numbers and they depend on the motor type and motor characteristics. They are: NV—output reduction ratio (ratio of the angular velocity at the output to the rotation speed at the input of the reductor), Nm—torque reduction ratio (ratio of the torques at the input and output of the reductor), JM—rotor torque, F—friction coefficient of the motor, Cn—motor’s mechanical constant (ratio of the motor shaft torque to the rotor current), RR—electrical resistance of the rotor circuit and CE—counter electromotive force coefficient. The electrical schema of the DC motors is represented in Figure 3. Based on (2) and (3) the mathematical model of the each actuator can be obtained. For the purpose of the further modeling it is necessary to express the equation of the actuator in the state space. We adopted the state vector as x̂=(θ,θ̇,IR)T. The order of the actuator’s equation is three. Without losing any precision the model of the actuator can be described with an order being equal to two because the inductance isLR≈0.  In this case, the nonlinear motor behavior is neglected.  A dynamic equation for the actuators is as follows: (4)x̂̇i=Aix̂i+biUi+fiMi,
where Aiis the motor’s system matrix having dimension 2×2, and biand fiare 2×1 vectors. MatricesAi, bi, and fican be represented as(5)Ai=[011-1JMi(FiNViNMi+CMiCEiRRi)],bi=[0CMiNViJMiRRi],fi=[0-1NViJMiNmi].Figure 2
5-DOF surgery module, absolute and local coordinate systems together with generalized coordinatesqi.Figure 3
Electric schema for DC motor. Each robot joint is motorized by direct current motor.The robotic system consists of the mechanical part (1) and actuators (4). The relationship between the actuator and manipulator’s mathematical model is twofold, via coordinates and via torques: θi=tiqi, wheretiis a transfer coefficient. For the described case, ti=1. Similarly, the actuator’s torque is equal to the torque on the jointMi=τi. Equation (1) can be rewritten as follows: (6)M(q)q̈+G(q,q̇)=τ,
where Mis the inertia matrix,Gis a matrix that represents the centrifugal, coriollis, and gravitational influence to the system, andτis the generalized torque vector. The state space vector of the whole system is adopted as(7)x=[x1x2⋯x2n-1x2n]T=[q1q̇1⋯qnq̇n]T,
wherenis the number of generalized coordinates, that is, number of DOF of the surgery module. Now, based on the equationθi=tiqi, it can be calculated for each actuator q̇i=Tix̂i and for the whole system q̇=Tx, keeping in mind that T=diag(Ti), Ti=(0ti). Based on (6) the mathematical equation of the mechanical part in the state space form is(8)τ=M(x)Tẋ+G(x).
Using the suggested relationships between the actuators and mechanical part via the coordinates and torques, the system model can be expressed by the following matrix equation:(9)ẋ=Ax+Bu+Fτ,
where A=diag(Ai) having order 2n, and B and F are 2n×n matrices described by the equations B=diag(bi) and F=diag(fi). Vector u is the input vectoru=(U1⋯Un)T.  Combining (6) and (8) and solving it in terms ofτthe following expression is obtained:(10)τ=(In-M(x)TF)-1[M(x)T(Ax+Bu)+G(x)].
The final form of the state space nonlinear dynamic equation is calculated by putting (10) into (9): (11)ẋ=Â(x)+B̂(x)u,
where the matrices from (11) are represented by (12): (12)Â(x)=(A+F(In-M(x)TF)-1M(x)TA)x+F(In-M(x)TF)-1G(x),B̂(x)=B+F(In-M(x)TF)-1M(x)TB.


### 3.2. System Dynamics for MPC

Assume that the contact surface is a scalar functionΦ:Rm→R1 (Figure 4):(13)Φ(p)=0.
The contact force vector is(14)f=DT(p)λ.λ is a scalar multiplier for the constraint function andD(p)is the gradient of the constraint function, described as in (15):(15)D(p)=∂Φ(p)∂p.
In the previous equation p∈ℜ3 is a position vector from the fixed reference frame to the constrained surface.  Additional constraint is p=Φ(q) at the contact point. The matrix function J(q) is defined as the Jacobian matrix function:(16)J(q)=∂Φ(q)∂q.
The influence of the contact force to the system can be described as(17)M(q)q̈+G(q,q̇)=τ+JT(q)f,
where M(q) denotes the inertia matrix function and G is the vector function which describes coriolis, centrifugal, and gravitational effects. q is a vector of the generalized coordinates and τ is the torque vector. Influence of the external contact forces is described by factorf. Combining (13)–(17) the mathematical model of the system is(18)M(q)q̈+G(q,q̇)=τ+JT(q)DT(M(q))λ(q,q̇,τ).
After an appropriate transformation, (18) is transformed to its matrix form:(19)[I000M(q)0000](q̇q̈λ̇)=[0I000JTDT000](qq̇λ)+[0-G(q,q̇)-M(q)Φ(p)]+[0I0]τ.
System (19) is linearized at the surrounding point where the needle is inserted into the patient.  The linearized mathematical model of the system is obtained, as follows:
(20)Eẋ(t)=Ax(t)+Bu(t)+d,
with singular matrix E∈ℜn×n, and vector d which represents the unmeasured disturbances such as needle deflection, tissue deformation, and displacement, A∈ℜn×n, B∈ℜm×r. x and u represent the system state-space vector and control vector, respectively, x∈ℜn, u∈ℜr.Figure 4
3D representation of the contact surface (prostate gland) which is variable and it is different for each case. Mathematical equation of contact surface is calculated using Euclidian robotic system software.The system defined by (20) is known as a singular system, descriptor system, semistate system, system of differential-algebraic equations, or generalized state space system. They arise naturally in many physical applications, such as electrical networks, aircraft and robotic systems, neutral delay and large-scale systems, economics, optimization problem, and constrained mechanics. The matrices of linearized system (20) are defined as follows:(21)A=[0I0∂∂q(G-JTDTλ)|00JTDT|0DJ|000],E=[I000M(q0)0000],B=[0I0],u=δτ,d=[0Δτ0].Now it is possible to find subspacesS and F, as in [34], having in mind that for state space X the system (20) can be represented as the direct sum of S and F,  X=S⊕F. As a result of the matrix transformation M applied to system (20), the slow and fast subsystems can be described by a mathematical formulation: (22)ẋs=Lsxs+Bsu+ds,Lfẋf=xf+Bfu+df,
with Ls=MA∣S, Lf=MA∣F, Bs=PMB, Bf=QMB, ds=PMd, and df=QMd for some linear transformations Q and P represented by matrices Q and P, respectively. For that case, the initial conditions are chosen to satisfy x0s=Px0 and x0f=Qx0. The solution of system (10) is x=xs+xfas in [35]: (23)xs=eLstx0s+∫0teLs(t-τ)Bsu(τ)dτ+∫0teLs(t-τ)d(τ)dτ,xf=-∑i=0v-1LfiBfui.Now it is possible to apply the NNPC and MPC controllers and to simulate dynamical behavior of the system when two types of the predictive approaches are considered.


## 3.1. System Dynamics for NNCP

Prior to implementing one of the suggested controllers to the system, it is essential to investigate the system behavior. Therefore, in this section, we have presented a mathematical model of the robotic system together with the mathematical model of the actuators. The derivation of dynamic equations is performed using the energy-based Lagrangian dynamic formulation. The mechanical part of the robotic subsystem is represented as a five-link open kinematic chain. All links are rigid bodies. We assumed that all materials for the system are homogenous. The generalized coordinates (q1, q2, q3, q4, q5) are adopted to describe each joint movement; see Figure 2. The general form of dynamic equation can be written as follows:
(1)∑α=1naαγ(q)q̈α+∑α=1n∑β=1nΓαβ,γ(q)q̇αq̇β=Qγg+Qγu,
where γ=1,…,5, aαγ(q) stands for the metrics tensor coefficients, Γαβ,γ(q) denotes the corresponding coefficients, and Qg, Qu are the generalized forces components. The former is a component due to gravitation and the latter is a component that corresponds to the actuator forces and torques. Each robot joint is motorized by a direct current motor, depending on the torque. Mathematical relationship between the output torque and the input motor torque for each motor is(2)NVNmJMθ̈+Fθ̇+M=CnNmIR.
Equation of the electrical equilibrium in the rotor circuit is: (3)LRİR+RRIR+CENVθ̇=U,
where θis the output rotation angle, Mis the output torque of the motor, IRis the rotor current, LRis the motor’s inductance, and U is the motor’s voltage. Other values are constant numbers and they depend on the motor type and motor characteristics. They are: NV—output reduction ratio (ratio of the angular velocity at the output to the rotation speed at the input of the reductor), Nm—torque reduction ratio (ratio of the torques at the input and output of the reductor), JM—rotor torque, F—friction coefficient of the motor, Cn—motor’s mechanical constant (ratio of the motor shaft torque to the rotor current), RR—electrical resistance of the rotor circuit and CE—counter electromotive force coefficient. The electrical schema of the DC motors is represented in Figure 3. Based on (2) and (3) the mathematical model of the each actuator can be obtained. For the purpose of the further modeling it is necessary to express the equation of the actuator in the state space. We adopted the state vector as x̂=(θ,θ̇,IR)T. The order of the actuator’s equation is three. Without losing any precision the model of the actuator can be described with an order being equal to two because the inductance isLR≈0.  In this case, the nonlinear motor behavior is neglected.  A dynamic equation for the actuators is as follows: (4)x̂̇i=Aix̂i+biUi+fiMi,
where Aiis the motor’s system matrix having dimension 2×2, and biand fiare 2×1 vectors. MatricesAi, bi, and fican be represented as(5)Ai=[011-1JMi(FiNViNMi+CMiCEiRRi)],bi=[0CMiNViJMiRRi],fi=[0-1NViJMiNmi].Figure 2
5-DOF surgery module, absolute and local coordinate systems together with generalized coordinatesqi.Figure 3
Electric schema for DC motor. Each robot joint is motorized by direct current motor.The robotic system consists of the mechanical part (1) and actuators (4). The relationship between the actuator and manipulator’s mathematical model is twofold, via coordinates and via torques: θi=tiqi, wheretiis a transfer coefficient. For the described case, ti=1. Similarly, the actuator’s torque is equal to the torque on the jointMi=τi. Equation (1) can be rewritten as follows: (6)M(q)q̈+G(q,q̇)=τ,
where Mis the inertia matrix,Gis a matrix that represents the centrifugal, coriollis, and gravitational influence to the system, andτis the generalized torque vector. The state space vector of the whole system is adopted as(7)x=[x1x2⋯x2n-1x2n]T=[q1q̇1⋯qnq̇n]T,
wherenis the number of generalized coordinates, that is, number of DOF of the surgery module. Now, based on the equationθi=tiqi, it can be calculated for each actuator q̇i=Tix̂i and for the whole system q̇=Tx, keeping in mind that T=diag(Ti), Ti=(0ti). Based on (6) the mathematical equation of the mechanical part in the state space form is(8)τ=M(x)Tẋ+G(x).
Using the suggested relationships between the actuators and mechanical part via the coordinates and torques, the system model can be expressed by the following matrix equation:(9)ẋ=Ax+Bu+Fτ,
where A=diag(Ai) having order 2n, and B and F are 2n×n matrices described by the equations B=diag(bi) and F=diag(fi). Vector u is the input vectoru=(U1⋯Un)T.  Combining (6) and (8) and solving it in terms ofτthe following expression is obtained:(10)τ=(In-M(x)TF)-1[M(x)T(Ax+Bu)+G(x)].
The final form of the state space nonlinear dynamic equation is calculated by putting (10) into (9): (11)ẋ=Â(x)+B̂(x)u,
where the matrices from (11) are represented by (12): (12)Â(x)=(A+F(In-M(x)TF)-1M(x)TA)x+F(In-M(x)TF)-1G(x),B̂(x)=B+F(In-M(x)TF)-1M(x)TB.


## 3.2. System Dynamics for MPC

Assume that the contact surface is a scalar functionΦ:Rm→R1 (Figure 4):(13)Φ(p)=0.
The contact force vector is(14)f=DT(p)λ.λ is a scalar multiplier for the constraint function andD(p)is the gradient of the constraint function, described as in (15):(15)D(p)=∂Φ(p)∂p.
In the previous equation p∈ℜ3 is a position vector from the fixed reference frame to the constrained surface.  Additional constraint is p=Φ(q) at the contact point. The matrix function J(q) is defined as the Jacobian matrix function:(16)J(q)=∂Φ(q)∂q.
The influence of the contact force to the system can be described as(17)M(q)q̈+G(q,q̇)=τ+JT(q)f,
where M(q) denotes the inertia matrix function and G is the vector function which describes coriolis, centrifugal, and gravitational effects. q is a vector of the generalized coordinates and τ is the torque vector. Influence of the external contact forces is described by factorf. Combining (13)–(17) the mathematical model of the system is(18)M(q)q̈+G(q,q̇)=τ+JT(q)DT(M(q))λ(q,q̇,τ).
After an appropriate transformation, (18) is transformed to its matrix form:(19)[I000M(q)0000](q̇q̈λ̇)=[0I000JTDT000](qq̇λ)+[0-G(q,q̇)-M(q)Φ(p)]+[0I0]τ.
System (19) is linearized at the surrounding point where the needle is inserted into the patient.  The linearized mathematical model of the system is obtained, as follows:
(20)Eẋ(t)=Ax(t)+Bu(t)+d,
with singular matrix E∈ℜn×n, and vector d which represents the unmeasured disturbances such as needle deflection, tissue deformation, and displacement, A∈ℜn×n, B∈ℜm×r. x and u represent the system state-space vector and control vector, respectively, x∈ℜn, u∈ℜr.Figure 4
3D representation of the contact surface (prostate gland) which is variable and it is different for each case. Mathematical equation of contact surface is calculated using Euclidian robotic system software.The system defined by (20) is known as a singular system, descriptor system, semistate system, system of differential-algebraic equations, or generalized state space system. They arise naturally in many physical applications, such as electrical networks, aircraft and robotic systems, neutral delay and large-scale systems, economics, optimization problem, and constrained mechanics. The matrices of linearized system (20) are defined as follows:(21)A=[0I0∂∂q(G-JTDTλ)|00JTDT|0DJ|000],E=[I000M(q0)0000],B=[0I0],u=δτ,d=[0Δτ0].Now it is possible to find subspacesS and F, as in [34], having in mind that for state space X the system (20) can be represented as the direct sum of S and F,  X=S⊕F. As a result of the matrix transformation M applied to system (20), the slow and fast subsystems can be described by a mathematical formulation: (22)ẋs=Lsxs+Bsu+ds,Lfẋf=xf+Bfu+df,
with Ls=MA∣S, Lf=MA∣F, Bs=PMB, Bf=QMB, ds=PMd, and df=QMd for some linear transformations Q and P represented by matrices Q and P, respectively. For that case, the initial conditions are chosen to satisfy x0s=Px0 and x0f=Qx0. The solution of system (10) is x=xs+xfas in [35]: (23)xs=eLstx0s+∫0teLs(t-τ)Bsu(τ)dτ+∫0teLs(t-τ)d(τ)dτ,xf=-∑i=0v-1LfiBfui.Now it is possible to apply the NNPC and MPC controllers and to simulate dynamical behavior of the system when two types of the predictive approaches are considered.


## 4. Results

### 4.1. Optimization Goals and Experiments

In the NNPC approach the ANN predicts insertion force. The adaptive optimization module was responsible to minimize the insertion force by recomputing the optimization parameters.In the MPC approach the slow subsystem was stabilized with state feedback(24)u=Kslxs.
A stabilizing feedback Ksl for the slow subsystem was obtained by minimization of the following expression:(25)min∫0∞(xxTQxs+uTPu)dt.
For data collection, we used the fully automated Robotic Brachytherapy System—EUCLIDIAN (Endo-Uro Computed Lattice for Intratumoral Delivery, Implantation, Ablation with Nanosensing); see Figure 1. The whole equipment was integrated into the system.  For cannula, we used Rev. A4, 304 Stainless steel needle (r = 2.98 mm ±.01 mm with 15° bevel angle) and for stylet  Rev. A4, 304 Stainless steel needle Vita Needle Company. The phantom was made of polyvinylchloride (PVC) plastic and PVC+ hardener in the ratio 80% to 20%, respectively (MF Manufacturing, TX). Force measurements were performed using two single-axis force sensors (Model 13, Honeywell Sensotech, Columbus, OH), one at each of the proximal ends of the stylet and cannula, and one 6-axis force-torque sensor (Nano17, ATI Industrial Automation, Apex, NC) at the distal end of the cannula. The deposited seeds were rounded stainless steel dummy seeds (BARD). The position of the needle tip and, consequently, the depth of deposition into the phantom were measured using optical encoders FAULHABER MicroMo HEDM5500J series 5500, attached to the cannula and stylet motors.


### 4.2. Force Prediction Using NNPC

Neural network controller uses a neural network model to predict future system responses to potential control signals (Figure5). The first task of the system was to train the neural network to represent the forward dynamics of the robotic system. An optimization algorithm then computes the control signals that optimize the future system parameters. For the MPC the dynamic model of the robotic system (11) is used to predict the future system behavior. An optimization algorithm is used to select the control input that optimizes the future performances. The error between the predicted force value and measured value is used as an ANN training signal.Figure 5
Predictive control system scheme: Neural network model is used to predict future plant performances. The controller calculates the control input that optimizes robotic system performances. After system identification procedure, robot model is used to predict future performances.A multilayer backpropagation neural network is employed to the NNPC scheme to model the nonlinear relationships between the robotic system input (actuators voltage) and the insertion force, in order to adapt the system to a variety of operating conditions and to acquire a more flexible learning ability. The differences of the predicted force and real-time measured values are within the margin of the 1.8% which may be satisfactory for the presented application, as in Figure6. The prediction error was shown in Figure 7. The time delay for such a complex nonlinear system could be modified if stronger system performaces are required. The predictive neural network had 3 hidden layers. The number of the time steps within the prediction errors that were minimized was 8. The number of the time steps within the control increments that are minimized was 2. In the presented simulation,  350 training epochs were used together with trainlm training function. It was found that the optimal update interval of 0.3s was good enough in the prediction procces, but our anticipation was that the suggested system performances might benefit if the future algorithms are applied in order to get a faster update rate.Figure 6
Insertion force prediction results using ANN for representative insertion: measured force (blue) and predicted force (green).Figure 7
Insertion force prediction error.In Figure8 the insertion force for both cases was represented: nonoptimized insertion and insertion when the ANN network approach is employed. It can be concluded that for optimized case the insertion force is smooth without peaks. Using the presented method, it is not possible to completely decrease the force during needle trajectory but high force gradient can be avoided. Consequently, it is possible that the tissue deformation and damage can be decreased.Figure 8
Optimized (green) and nonoptimized (blue) insertion forces measured by force sensor installed in the robotic system.


### 4.3. Insertion Force Control and MPC

The purpose of this approach is to control the reactive insertion force on one side and to predict and compensate the impact of the measured and unmeasured disturbances rather than waiting until the effect appears at the output of the system on the other side. Because of the fact that this procedure required an accurate robot system model, it was necessary to obtain a more precise model. That is a reason for adopting the robotic system model as a singular system of differential equations. In [33] it is shown that the impulsive behavior of the system can be avoided using appropriate initial conditions, defined by x0s=Px0and x0f=Qx0. By using the described approach the fast subsystem will not induce impulsive behavior. Moreover, it can be concluded as stated previously and from (23) that there is little need to find fast feedback to eliminate the impulsive behavior. The necessary task was to find an appropriate feedback for the slow subsystem. The control scheme for the slow subsystem is represented in Figure 9, as suggested in [34].Figure 9
Insertion force control scheme for the slow subsystem. Disturbanced is given to MPC and its effect is predicted and compensated before the effect appears at the system output.Furthermore, when the MPC controller is applied, the main objective is to hold the insertion force at the predefined acceptable reference value, by adjusting the control signal from actuators in order to minimize the prostate displacement, that is, the reactive force which acts upon the tissue. Using this approach it is possible to decrease the insertion force during insertion trajectory. The needle displacement is an unmeasured disturbance and the controller provides feedback compensation for such disturbances. For the insertion force the controller provides feedforward compensation. Various noise effects can corrupt the measurements. The noise could vary randomly with a zero mean or could exhibit a nonzero, drifting bias. The MPC uses a filtering method for removing estimated noise component. At the beginning of each sampling instant, the controller estimates the current system state. Accurate knowledge of the state improves prediction accuracy which improves controller performances. As a result of passive insertion force control, multiplierλchanges values during the insertion and the force is minimized based on (25). During the insertion passive force control becomes active and keeps passive insertion force close to the predefined minimal value. When the MPC control approach was implemented, it is possible to decrease the insertion force, as it is shown in Figure 10. Also, peaks during the insertion are minimized. These conclusions give us a reason to believe that tissue deformation can be decreased better than using the NNPC approach. The previous hypothesis was proven after the experiments, as in Figure 11. Figure 11 shows the result of tissue deformation for both insertion approaches, the MPC and classic insertion with open loop control. It can be observed that the open-loop controller is not able to predict insertion force and consequently to minimize tissue displacement. A system with the MPC has better optimization characteristics and faster response. Consequently, tissue deformation is less compared to the open-loop case as well as to the NNPC approach. For the needle displacement, the results are similar. As a result of the experiments, the needle displacement, and therefore seed deposition precision, for insertion depth of 8 cm in the first case (NNPC) was between 0.8 and 1.1 mm. For the MPC it was 0.2–0.4 mm. The results in the second approach were better because the MPC includes an influence on the needle deflection as an unmeasured disturbance. Presented feedback compensation method resulted in needle deflection which ensures clinically acceptable precision for radioactive seed deposition (less then ±2 mm, according to [36]).Figure 10
Controlled force (red) and insertion force when MPC is not implemented (blue).Tissue displacements for the MPC (green) and for open-loop control approach (red). Disturbanced is given to MPC and its effect is predicted and compensated before effect appears at the system output.
(a)(b)(c)(d)


## 4.1. Optimization Goals and Experiments

In the NNPC approach the ANN predicts insertion force. The adaptive optimization module was responsible to minimize the insertion force by recomputing the optimization parameters.In the MPC approach the slow subsystem was stabilized with state feedback(24)u=Kslxs.
A stabilizing feedback Ksl for the slow subsystem was obtained by minimization of the following expression:(25)min∫0∞(xxTQxs+uTPu)dt.
For data collection, we used the fully automated Robotic Brachytherapy System—EUCLIDIAN (Endo-Uro Computed Lattice for Intratumoral Delivery, Implantation, Ablation with Nanosensing); see Figure 1. The whole equipment was integrated into the system.  For cannula, we used Rev. A4, 304 Stainless steel needle (r = 2.98 mm ±.01 mm with 15° bevel angle) and for stylet  Rev. A4, 304 Stainless steel needle Vita Needle Company. The phantom was made of polyvinylchloride (PVC) plastic and PVC+ hardener in the ratio 80% to 20%, respectively (MF Manufacturing, TX). Force measurements were performed using two single-axis force sensors (Model 13, Honeywell Sensotech, Columbus, OH), one at each of the proximal ends of the stylet and cannula, and one 6-axis force-torque sensor (Nano17, ATI Industrial Automation, Apex, NC) at the distal end of the cannula. The deposited seeds were rounded stainless steel dummy seeds (BARD). The position of the needle tip and, consequently, the depth of deposition into the phantom were measured using optical encoders FAULHABER MicroMo HEDM5500J series 5500, attached to the cannula and stylet motors.


## 4.2. Force Prediction Using NNPC

Neural network controller uses a neural network model to predict future system responses to potential control signals (Figure5). The first task of the system was to train the neural network to represent the forward dynamics of the robotic system. An optimization algorithm then computes the control signals that optimize the future system parameters. For the MPC the dynamic model of the robotic system (11) is used to predict the future system behavior. An optimization algorithm is used to select the control input that optimizes the future performances. The error between the predicted force value and measured value is used as an ANN training signal.Figure 5
Predictive control system scheme: Neural network model is used to predict future plant performances. The controller calculates the control input that optimizes robotic system performances. After system identification procedure, robot model is used to predict future performances.A multilayer backpropagation neural network is employed to the NNPC scheme to model the nonlinear relationships between the robotic system input (actuators voltage) and the insertion force, in order to adapt the system to a variety of operating conditions and to acquire a more flexible learning ability. The differences of the predicted force and real-time measured values are within the margin of the 1.8% which may be satisfactory for the presented application, as in Figure6. The prediction error was shown in Figure 7. The time delay for such a complex nonlinear system could be modified if stronger system performaces are required. The predictive neural network had 3 hidden layers. The number of the time steps within the prediction errors that were minimized was 8. The number of the time steps within the control increments that are minimized was 2. In the presented simulation,  350 training epochs were used together with trainlm training function. It was found that the optimal update interval of 0.3s was good enough in the prediction procces, but our anticipation was that the suggested system performances might benefit if the future algorithms are applied in order to get a faster update rate.Figure 6
Insertion force prediction results using ANN for representative insertion: measured force (blue) and predicted force (green).Figure 7
Insertion force prediction error.In Figure8 the insertion force for both cases was represented: nonoptimized insertion and insertion when the ANN network approach is employed. It can be concluded that for optimized case the insertion force is smooth without peaks. Using the presented method, it is not possible to completely decrease the force during needle trajectory but high force gradient can be avoided. Consequently, it is possible that the tissue deformation and damage can be decreased.Figure 8
Optimized (green) and nonoptimized (blue) insertion forces measured by force sensor installed in the robotic system.


## 4.3. Insertion Force Control and MPC

The purpose of this approach is to control the reactive insertion force on one side and to predict and compensate the impact of the measured and unmeasured disturbances rather than waiting until the effect appears at the output of the system on the other side. Because of the fact that this procedure required an accurate robot system model, it was necessary to obtain a more precise model. That is a reason for adopting the robotic system model as a singular system of differential equations. In [33] it is shown that the impulsive behavior of the system can be avoided using appropriate initial conditions, defined by x0s=Px0and x0f=Qx0. By using the described approach the fast subsystem will not induce impulsive behavior. Moreover, it can be concluded as stated previously and from (23) that there is little need to find fast feedback to eliminate the impulsive behavior. The necessary task was to find an appropriate feedback for the slow subsystem. The control scheme for the slow subsystem is represented in Figure 9, as suggested in [34].Figure 9
Insertion force control scheme for the slow subsystem. Disturbanced is given to MPC and its effect is predicted and compensated before the effect appears at the system output.Furthermore, when the MPC controller is applied, the main objective is to hold the insertion force at the predefined acceptable reference value, by adjusting the control signal from actuators in order to minimize the prostate displacement, that is, the reactive force which acts upon the tissue. Using this approach it is possible to decrease the insertion force during insertion trajectory. The needle displacement is an unmeasured disturbance and the controller provides feedback compensation for such disturbances. For the insertion force the controller provides feedforward compensation. Various noise effects can corrupt the measurements. The noise could vary randomly with a zero mean or could exhibit a nonzero, drifting bias. The MPC uses a filtering method for removing estimated noise component. At the beginning of each sampling instant, the controller estimates the current system state. Accurate knowledge of the state improves prediction accuracy which improves controller performances. As a result of passive insertion force control, multiplierλchanges values during the insertion and the force is minimized based on (25). During the insertion passive force control becomes active and keeps passive insertion force close to the predefined minimal value. When the MPC control approach was implemented, it is possible to decrease the insertion force, as it is shown in Figure 10. Also, peaks during the insertion are minimized. These conclusions give us a reason to believe that tissue deformation can be decreased better than using the NNPC approach. The previous hypothesis was proven after the experiments, as in Figure 11. Figure 11 shows the result of tissue deformation for both insertion approaches, the MPC and classic insertion with open loop control. It can be observed that the open-loop controller is not able to predict insertion force and consequently to minimize tissue displacement. A system with the MPC has better optimization characteristics and faster response. Consequently, tissue deformation is less compared to the open-loop case as well as to the NNPC approach. For the needle displacement, the results are similar. As a result of the experiments, the needle displacement, and therefore seed deposition precision, for insertion depth of 8 cm in the first case (NNPC) was between 0.8 and 1.1 mm. For the MPC it was 0.2–0.4 mm. The results in the second approach were better because the MPC includes an influence on the needle deflection as an unmeasured disturbance. Presented feedback compensation method resulted in needle deflection which ensures clinically acceptable precision for radioactive seed deposition (less then ±2 mm, according to [36]).Figure 10
Controlled force (red) and insertion force when MPC is not implemented (blue).Tissue displacements for the MPC (green) and for open-loop control approach (red). Disturbanced is given to MPC and its effect is predicted and compensated before effect appears at the system output.
(a)(b)(c)(d)


## 5. Conclusion and Future Work

In this article, we used the neural network predictive control (NNPC) to predict an insertion force. The NNPC uses the nonlinear model of the robotic system to predict future system performances. The controller then calculates the control input vector (motors voltages) that optimizes the robot performances based on predicted force values over finite time interval. The control system is designed to maintain practical stability even when the force gradient is high. It is shown that the neural network approach is suitable for improving performances of the robotic system with complex dynamic behavior.In the second approach the MPC method was applied. This method shows better results than the NNPC method. There are two reasons for that. First is that the mathematical model of the system used for NNPC does not take the insertion force into account. The insertion force was measured but the robot dynamics does not depend directly on the insertion force. Based on the experiments, the results show that this method has a slower response, and sometimes it is not good enough for displacement compensation. Actually, the compensation of the displacement for the NNCP method is performed by decreasing force gradient. It was a consequence of the insertion force prediction and optimization of the insertion procedure. That is a second reason why MPC was more suitable for this kind of application. On the other hand, the mathematical model in the MPC case includes the insertion force. Moreover, the system is described by using a singular system of differential equations. That approach gives a more accurate mathematical model which is one of the basic prerequisites for the MPC method. An influence of the insertion force can be changed by observing changes of the Lagrange multiplierλ. Here, not only the control of the passive insertion force (the reaction force) were implemented, but also the disturbances are taken into account. All disturbances were divided into two groups, measured and unmeasured. The insertion force was a measured disturbance to the system and the needle and tissue displacement belonged to another group of disturbances. Together with the control of the passive force, position control and velocity control were applied as well.The second approach (MPC) can include more procedure specific criteria such as needle acceleration or rotation for each specific encounter, or patent specific criteria such as age, prostate volume, activity, or even prostate density or texture. Based on the presented results an effective motion planning algorithm will be developed as well as the improved predictive and optimization algorithm.



---

*Source: 581840-2010-05-20.xml*