# The Design of Compact Robotic-Assisted Needle Position System with MPC-Based Remote Control

Jing Guo, Yi Liu, Jin Wang, Chao Zeng, Jie Huang, Chao Liu

► To cite this version:

Jing Guo, Yi Liu, Jin Wang, Chao Zeng, Jie Huang, et al. The Design of Compact Robotic-Assisted Needle Position System with MPC-Based Remote Control. Complexity, 2020, 2020, pp.#5406084. 10.1155/2020/5406084. lirmm-03287026

HAL Id: lirmm-03287026

https://hal-lirmm.ccsd.cnrs.fr/lirmm-03287026

Submitted on 15 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Article

# The Design of Compact Robotic-Assisted Needle Position System with MPC-Based Remote Control

Jing Guo,<fnref n="1" /> Yi Liu,<fnref n="1" /> Jin Wang,<fnref n="1" /> Chao Zeng ID,<fnref n="1" /> Jie Huang ID,<fnref n="2" /> and Chao Liu<fnref n="3" />

Correspondence should be addressed to Chao Zeng; auczeng@gdut.edu.cn and Jie Huang; jieh2015@yahoo.com

Received 28 June 2020; Revised 7 August 2020; Accepted 30 August 2020; Published 16 September 2020

Academic Editor: Yan Wu

Copyright © 2020 Jing Guo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article introduces the design and control performance of a lightweight, flexible, 4-degree-of-freedom (DOF) parallel robot for percutaneous biopsy guided by computed tomography (CT). At present, the CT guidance method allows surgeons to quickly locate the lesion area; however, it is necessary to manually adjust the position of the puncture needle for insertion. In this paper, a three-dimensional assisted method is used to infer the control input required to reach the target point through the kinematic model of the robot. A Kalman filter is designed to estimate model parameters and obtain a more accurate model. To further improve the control performance of the robot system, a model-based control method—the model predictive control (MPC) controller—is used to increase the accuracy of the needle position in the developed robot system. In this way, medical efficiency is improved while reducing the burden on the surgeon.

## 1. Introduction

The surgeries have been advanced in the past decades thanks to novel surgical techniques [1], such as laparoscopy surgeries, in which a surgeon can operate a surgery through tiny holes on a human body for removing pathological tissue. However, the hand-hold tools still make the operation hard to be performed due to limited ergonomics and accuracy. Meanwhile, the surgeons may feel fatigued when operating for long time during the surgery. In order to improve the quality of surgical operation, the robot assistant technique has been introduced into the operation room which has eased the burden of the surgeon during the surgery [2-4]. Nowadays, the robotic systems have been developed and used for various surgeries such as biopsy, brachytherapy, and tumor ablation [5-7].

Before performing a minimally invasive surgical (MIS) task, the coordinates of the robot are usually registered to the medical image coordinate system, so that the robot is controlled by analysis of the lesion's position for further

operations. According to the characteristics of MIS, the pose of the robot end-effector before needle insertion during surgery becomes very important. The needle positioning system of the surgical assistant robot can be independent of the insertion action [8-10]. The surgeon can guide the robot's needle to the insertion position of the lesion on the basis of the medical image through image guidance [11-16] and then manually insert the needle into the lesion area. Medical imaging, acting as a real-time feedback tool for the needle pose, is crucial in the process of needle positioning and control. In current needle-guided systems, CT [11-13, 17, 18] is often used. Today's commercially available minimally invasive robotic systems (MIRS) maintain and expand the flexibility of the surgeon's hand [19-22]. Navarro-Alarcon et al. [12] developed a 3-DOF needle driver for biopsy. The positioning and insertion of the needle are achieved through three interfaces, the first two of which are aligned to the target, and the third is inserted. Koseki et al. [23] established a cooperative manipulation structure which uses the optical linear encoder to measure the needle

<footnote n="1">¹School of Automation, Guangdong University of Technology, Guangzhou, China</footnote>

<footnote n="2">²Guangdong Lung Cancer Institute, Guangdong Provincial Key Laboratory of Translational Medicine in Lung Cancer, Guangdong Provincial People's Hospital, Guangdong Academy of Medical Sciences, Guangzhou, China</footnote>

<footnote n="3">³Department of Robotics, LIRMM-CNRS, Montpellier, France</footnote>
position. However, manual remote operation control requires multiple scans using CT to obtain accurate targets. The robot-assisted positioning approach can locally record the workspace and medical images of the robot system to reduce the risk of cancer induced by repeated CT scans [24]. This approach poses challenges to the compatibility and size of the robot system.

In the development of a surgical robot, another major problem is the network system in the remote operating system, which is caused by long distance or wireless link [23–27]. Sampling and delay are involved in this problem, which is a challenge for remote controlling of the robot. The solution to the delay problem in remote is discussed in [28–31]. However, the local optimizer on the slave side of the robot will also cause a delay [32, 33]. This part of the delay will also affect the remote operating system following the feedback signal. Dong et al. [34] modified standard functions to improve the response speed. Norizuki and Uchimura [35] embedded the MPC controller into the remote operating system to reduce the impact of the optimizer's delay in the system.

In order to improve the quality of the surgery, the accuracy of the needle pose becomes crucial. The needle pose is affected by the robot control system and medical images. However, factors such as friction in the design of the robot structure and noise in the feedback of medical images are inevitable, which makes it difficult to accurately control the needle positioning [36, 37]. Thus, the closed-loop control method becomes a good solution [38, 39]. In this paper, firstly, as shown in Figure 1, a novel light-weighted puncture robotic system is provided. A surgeon controls the medical robot through the master in the remote system. The medical robot can provide feedback of the end-effector pose and force to the surgeon. The size of the designed robot is 16 cm high, 37.6 cm long, and 37.6 cm wide. In addition, the robot is controlled by four motors, which leads to the coupling of robot control. Due to the characteristics of the robot structure, specific constraints are proposed to decouple the kinematics. Finally, we develop a related needle positioning system on this robotic platform, which uses an MPC strategy to achieve the robotic local closed-loop control of the needle tip position. This approach not only solves the problem of the time delay (TD) of the needle positioning system in local, but also reduces the cost and time of personnel training required for the operation of the robotic needle positioning system.

The structure of the paper is organized as follows: in Section 2, the overall robotic system design is introduced, including the derivation of constrained kinematic decoupling approach. Section 3 introduces model parameter estimation. Section 4 describes MPC control method. In Section 5, the experiments and the analysis of the experimental results are introduced. We conclude this article in Section 6.

In our previous work [40], we have designed a remote system which contains a master and a slave robot (see Figure 1). A surgeon can control the slave robot through the master side, and during a task, the desired position of the end-effector can be marked and recorded. The end-effector

<figure><img src="image_2.png" /><figcaption>FIGURE 1: Surgical assistant robot remote operating system.</figcaption></figure>

position of the slave robot is updated via medical imaging or other sensors, and such position feedback is sent to master. There are many ways to solve the remote communication scheme between the master and slave robots [28–31]. However, when the slave robot uses the optimized control method to improve the control effect, a delay will occur in the slave robot, which will affect the use of the remote operating system. Thus, the local controller design of the slave robot needs to be considered.

Because the robotic control system is a multi-input multioutput (MIMO) system, the controller selected for the slave robot should be able to adapt to the MIMO system. The proposed slave robot control system can be seen in Figure 2(a). The surgeon uses CT to mark the target lesion area or trajectory in the robot workspace. The surgical robot obtains the control trajectory which serves as the reference input for the MPC controller through kinematic model. The MPC controller optimizes the current control input to the surgical robot according to the reference trajectory and the reference state of the surgical robot. Furthermore, in order to improve the control performance, the state-space model of the robot required by the MPC controller should be accurate. In this paper, according to the structure of the surgical robot, the kinematics model is divided into a linear part and a nonlinear part. The linear part uses the model parameter estimator to identify the model parameters, and the nonlinear part is linearized using the Taylor formula. As shown in Figure 2(b), the model parameter estimator of the linear part is composed of an order judger and a parameter estimator. The order judger is used to judge the order of the linear system. The Kalman filter is selected as the parameter estimator to estimate the parameters of the linear system.

## 2. Kinematic Modeling

In this section, the modeling process of surgical robot kinematics will be described. The robotic kinematic is analyzed for the establishment of model parameter estimator and MPC controller. A suitable kinematics model is constructed by analyzing the mechanical structure of the parallel puncture needle robot [40]. Figure 3 shows the three-dimensional mechanical structure and sectional structure of the surgical robot. The blocks, i.e., block-a and block-b, can be moved by controlling the angle of four motors. When the
<figure><img src="image_3.png" /><figcaption>FIGURE 2: (a) The overall control framework of the surgical robot positioning system. (b) A schematic of the model parameter estimator of the linear part of the kinematics model.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3: The three-dimensional mechanical structure and sectional structure of the surgical robot. It uses two parallel blocks to control the insertion direction and position of the same puncture needle.</figcaption></figure>

positions of the two blocks change, the position and insertion direction of the puncture needle will change accordingly. Among this structure, the vertical distance *d* between the two blocks and the length *h* of the puncture needle below the block-*b* are both constant.

In order to obtain a suitable kinematics model, the coordinate relationship as shown in Figure 4 is established. The fixed world coordinate system is at the initial center position of block-*b*. $^wT_a$ is the coordinate system of block-*a*, which is used to describe the moving distance of block-*a*. Similarly, the value of the coordinate system $^wT_b$ is the moving distance of block-*b*. Coordinate system $^wT_e$ indicates the position of the tip of the puncture needle.

**2.1. Forward Kinematic.** The purpose of forward kinematics is to obtain the position of the tip of the puncture needle $x_e$, $y_e$, and $z_e$, and the forward kinematics of the surgical robot can be obtained using geometric theorem. As shown in Figure 4, firstly, the rotation of the motor will drive the block. The moving distance of the two blocks $x_a$, $x_b$, $y_a$, and $y_b$ will change the angle of the puncture needle $\theta$.

<figure><img src="image_5.png" /><figcaption>FIGURE 4: The surgical robot's coordinate relationship of the kinematic model. Define the world coordinate system and the coordinate system of each joint to describe the kinematics model.</figcaption></figure>

$$x_a = \rho \beta_{a1}, \qquad (1)$$

$$\theta = \arctan\left(\frac{\sqrt{(y_a - y_b)^2 + (x_a - x_b)^2}}{d}\right), \qquad (2)$$

$$\alpha = \arctan\left(\frac{x_a - x_b}{y_a - y_b}\right), \qquad (3)$$

where $\rho$ is the pitch of the screw connected to the block; $\beta$ is the rotation angle of the motor, and the relationship between the rotation angle of other motors and its block coordinates
is similar; $\alpha$ is the angle generated by two blocks on the X-Y plane; next, for coordinate system ${}^wT_b$}, the change in the position of the tip of the puncture needle $\Delta x$, $\Delta y$, and $\Delta z$ can be described by $\theta$, $\alpha$, and the constant length of the needle $h$:

$$ \Delta x = h \cdot \sin(\theta) \cdot \sin(\alpha), \quad (4) $$

$$ \Delta y = h \cdot \sin(\theta) \cdot \cos(\alpha), \quad (5) $$

$$ \Delta z = h \cdot \cos(\theta). \quad (6) $$

Finally, the position of the puncture needle's tip $x_e$, $y_e$, and $z_e$ relative to the fixed world coordinate system can be obtained as follows:

$$ x_e = \Delta x + x_b, \quad (7) $$

$$ y_e = \Delta y + y_b, \quad (8) $$

$$ z_e = \Delta z. \quad (9) $$

**2.2. Inverse Kinematic.** The surgical robot is a parallel structure; hence, there is a coupling relationship between the control inputs. To address this problem, constraints are first defined to decouple the kinematic model. From the point of view of the positioning range of the puncture needle, it is better to have a larger range that the puncture needle can reach in a limited space. Thus, the first constraint is that the two blocks should move in opposite directions. On the other hand, from the view of the lesion area marked by surgeon, it can be seen that the position of the tip of the puncture needle in the X-Y plane is of more significance. In order to obtain a feasible solution for the motor control input, the second constraint is defined. That is, when block-a moves to the limit position and still cannot reach the lesion area, block-b is then moved. The inverse kinematics model is constructed using geometric theorem based on the constraints:

$$ \begin{gathered} [x_b, y_b] = [0, 0] \\ \text{s.t.} \quad -\eta \le x_t, y_t \le \eta, \end{gathered} \quad (10) $$

and

$$ \begin{gathered} [x_a, y_a] = [-2\operatorname{sgn}(x_t), -2\operatorname{sgn}(y_t)] \\ \text{s.t.} \quad |x_t, y_t| \le -\eta, \end{gathered} \quad (11) $$

where $\eta$ is the largest distance that the tip of the puncture needle can reach when only block-a is in motion; $x_t$ and $y_t$ are the lesion's target point coordinates marked by the surgeon. However, no complete feasible solution can be obtained yet. Then, the solution of the coordinates of block-a within the limit $\eta$ is discussed:

$$ \begin{bmatrix} x_a \\ y_a \end{bmatrix} = \begin{bmatrix} -d \cdot \operatorname{sgn}(x_t) \tan(\arcsin(\theta)) \sin\left(\arctan\left(\frac{x_t}{y_t}\right)\right) \\ -d \cdot \operatorname{sgn}(y_t) \tan(\arcsin(\theta)) \cos\left(\arctan\left(\frac{x_t}{y_t}\right)\right) \end{bmatrix}, \quad (12) $$

where $\theta$ is the angle between the block and puncture needle:

$$ \theta = \arccos\left(\frac{-z_t}{h}\right). \quad (13) $$

So far, a feasible solution within the limit $\eta$ from equations (10) and (12) is obtained:

$$ \begin{bmatrix} x_a \\ y_a \\ x_b \\ y_b \end{bmatrix} = \begin{bmatrix} -d \cdot \operatorname{sgn}(x_t) \tan(\arcsin(\theta)) \sin\left(\arctan\left(\frac{x_t}{y_t}\right)\right) \\ -d \cdot \operatorname{sgn}(y_t) \tan(\arcsin(\theta)) \cos\left(\arctan\left(\frac{x_t}{y_t}\right)\right) \\ 0 \\ 0 \end{bmatrix} $$ 
 s.t. $-\eta \le x_t, y_t \le \eta$. \quad (14)

Finally, the feasible solution in the other case will be analyzed to obtain the complete feasible solution. When the lesion area exceeds the limit $\eta$, the coordinate values of block-a are fixed. Therefore, according to the angle $\theta$ between the block and the puncture needle, the relationship between $x_b$ and $x_b$ can be acquired:

$$ \sqrt{(x_a - x_b)^2 + (y_a - y_b)^2} = d \cdot \tan(\theta), \quad (15) $$

$$ \frac{x_a - x_b}{y_a - y_b} = \frac{x_t - x_b}{y_t - y_b}. \quad (16) $$

The solution of coordinate value of block-b can be acquired by equations (15) and (16):

$$ \begin{bmatrix} x_b \\ y_b \end{bmatrix} = \begin{bmatrix} \frac{-EB - \sqrt{EB^2 - 4 \cdot EA \cdot EC}}{2 \cdot EA} \\ \frac{E - 2 \cdot (x_a - x_t)}{y_a - y_t} \end{bmatrix}, \quad (17) $$

where E, EA, EB, and EC are known parameters as follows:

$$ \begin{align} E &= \frac{h^2 - z_t^2 + x_a^2 - y_t^2 + y_a^2}{2}, \notag \\ EA &= \frac{(x_a - x_t)^2}{(y_a - y_t)^2} + 1, \notag \\ EB &= \frac{2d(y_a - y_t)(x_a - x_t) - 2(y_a - y_t)^2 x_t - 2E(x_a - x_t)}{(y_a - y_t)^2}, \notag \\ EC &= \frac{(y_a - y_t)^2 (x_t^2 + y_t^2 - h^2 + z_t^2) - 2dE(y_a - y_t)}{(y_a - y_t)^2}. \tag{18} \end{align} $$

Then, the feasible solution that exceeds the limit $\eta$ is obtained by equations (11) and (17):
$$ 
\begin{bmatrix} x_a \\ y_a \\ x_b \\ y_b \end{bmatrix} = \begin{bmatrix} -2\operatorname{sgn}(x_t) \\ -2\operatorname{sgn}(y_t) \\ \frac{-EB - \sqrt{EB^2 - 4 \cdot EA \cdot EC}}{2 \cdot EA} \\ \frac{E - 2 \cdot (x_a - x_t)}{y_a - y_t} \end{bmatrix} \quad (19)
$$

$$ \text{s.t.} \quad |x_t, y_t| \leq -\eta. $$

Now, the complete solution of the kinematic model of the surgical robot has been obtained by equations (14) and (19). The corresponding motor control angle can be obtained through the lesion area marked by the surgeon.

### 3. Model Parameter Estimation

After the desired control input is obtained, it cannot, in fact, directly act on the robotic controller. In order to make the state-space model of the robot closer to the actual model, the transmission process of the surgical robot will be divided into two parts for modeling. The first part is the linear system model of the block movement, and the second part is the nonlinear system model of the puncture needle movement. The determinant ratio is used to estimate the order of the

linear system model, and a Kalman filter is established to estimate the linear system model parameters. Taylor formula is used to linearize the nonlinear system model into the controller.

**3.1. Linear Part.** Firstly, the linear system model is a MIMO system; however, each input corresponds to only one output. Thus, the MIMO system is equivalent to four single-input single-output (SISO) systems. The order and parameters of the SISO system are considered:

$$ y(k) + a_1 y(k-1) + \cdots + a_n y(k-n) = b_0 u(k) + e(k), \quad (20) $$

where $e(k)$ is the Gaussian white noise with mean zero, $u(k)$ is the input signal, and $y(k)$ is the output signal. The input and output signals at each moment are known, but parameters are unknown. Then, input signal and output signal are used to construct determinants for comparison to obtain the ratio $DR^*(\hat{n})$:

$$ DR^*(\hat{n}) = \frac{\det[\mathbf{H}^*(\hat{n})]}{\det[\mathbf{H}^*(\hat{n}+1)]}, \quad (21) $$

where $\hat{n}$ is the estimated order and $\mathbf{H}^*$ is a matrix composed of input and output and has the following relationship:

$$ \mathbf{H}_{\hat{n}}^* = \begin{bmatrix} y(\hat{n}) & y(\hat{n}-1) & \cdots & y(1) & u(\hat{n}) & u(\hat{n}-1) & \cdots & u(1) \\ y(\hat{n}+1) & y(\hat{n}) & \cdots & y(2) & u(\hat{n}+1) & u(\hat{n}) & \cdots & u(2) \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ y(\hat{n}+L-1) & y(\hat{n}+L-2) & \cdots & y(L) & u(\hat{n}+L-1) & u(\hat{n}+L-2) & \cdots & u(L) \end{bmatrix}, \quad (22) $$

$$ \mathbf{H}^*(\hat{n}) = \frac{1}{L} \mathbf{H}_{\hat{n}}^* \mathbf{H}_{\hat{n}}^*, $$

where $L$ is the sampling length of the signal. $\hat{n}$ gradually increases from one. If only $DR^*(\hat{n})$ has a significant growth compared to $DR^*(\hat{n}-1)$, $\hat{n}$ is the estimated order.

Next, a Kalman filter is designed to dynamically identify the system's unknown parameters $a_i$ ($i = 1, 2, \dots, n$) and $b_0$. The Kalman filter is usually used to optimally estimate the state value of the state-space equation, which is composed of the unknown parameters:

$$ \begin{cases} \mathbf{X}(k+1) = \mathbf{\Phi}\mathbf{X}(k) + \Gamma\mathbf{W}(k), \\ \mathbf{Y}(k) = \mathbf{H}(k)\mathbf{X}(k) + \mathbf{V}(k), \end{cases} \quad (23) $$

where both $\mathbf{V}(k)$ and $\mathbf{W}(k)$ are noise matrix composed of Gaussian white noise with zero mean. $\Gamma$ is the process noise drive matrix. $\mathbf{X}(k)$ is the state matrix at the $k$-th sampling moment. $\mathbf{\Phi}$ is the state transition matrix at the $k$-th sampling moment. $\mathbf{Y}(k)$ is the output matrix at the $k$-th sampling moment. And we assume that the variances of the two noises

are $\mathbf{Q}$ and $\mathbf{R}$. Then, the update equation of Kalman filter can be derived based on Euclidean theorem:

$$ \left\{ \begin{aligned} \hat{\mathbf{X}}(k+1|k) &= \mathbf{\Phi}\hat{\mathbf{X}}(k|k), \\ \hat{\mathbf{X}}(k+1|k+1) &= \hat{\mathbf{X}}(k+1|k) + \mathbf{K}(k+1)(\mathbf{Y}(k+1) - \mathbf{H}\hat{\mathbf{X}}(k+1|k)), \\ \mathbf{K}(k+1) &= \mathbf{K}(k+1|k)\mathbf{H}^T[\mathbf{H}\mathbf{P}(k+1|k)\mathbf{H}^T + \mathbf{R}]^{-1}, \\ \mathbf{P}(k+1|k) &= \mathbf{\Phi}\mathbf{P}(k|k)\mathbf{\Phi}^T + \Gamma\mathbf{Q}\Gamma^T, \\ \mathbf{P}(k+1|k+1) &= [\mathbf{I}_n - \mathbf{K}(k+1)\mathbf{H}^T]\mathbf{P}(k+1|k), \end{aligned} \right. \quad (24) $$

where $\hat{\mathbf{X}}(k|k)$ is the optimal estimate at the $k$-th sampling moment, $\mathbf{P}(k|k)$ is the error covariance of the current sampling time, and $\mathbf{K}(k+1|k)$ is the gain of the Kalman filter, which is used to update the values of $\hat{\mathbf{X}}(k|k)$ and
$\mathbf{P}(k|k)$. In this paper, the estimated state matrix with unknown parameters is constructed:

$$ 
\left\{ 
\begin{array}{l}
x_1(k) = a_1(k), \\
x_2(k) = a_2(k), \\
\vdots \\
x_n(k) = a_n(k), \\
x_{n+1}(k) = b_0(k).
\end{array}
\right. 
\qquad (25) $$

Then the system's state equation is obtained as follows:

$$ \mathbf{X}(k+1) = \mathbf{X}(k) + \mathbf{W}(k), \qquad (26) $$

where $\mathbf{X}(k)$ is composed of $x_i(k)$ ($i = 1, 2, \dots, n+1$). We use the observed input and output data to construct the matrix $\mathbf{H}(k)$ and state equation:

$$ 
\begin{aligned}
\mathbf{H}(k) = & [-y(k-1), -y(k-2), \dots, -y(k-n), u(k)], \\
\begin{cases}
\mathbf{X}(k+1) = \mathbf{X}(k) + \mathbf{W}(k), \\
y(k) = \mathbf{H}(k)\mathbf{X}(k) + e(k).
\end{cases}
\end{aligned} 
\qquad (27) $$

Finally, state equation (28) can be updated by equation (25) to obtain the actual system parameters. The SISO systems of four motors are combined to obtain the MIMO system of the linear part:

$$ 
\mathbf{A} = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ -a_n & -a_{n-1} & -a_{n-2} & \cdots & -a_1 \end{bmatrix}, \qquad (28) $$

$$ \mathbf{B}^T = [0 \ 0 \ \cdots \ 0 \ b_0], \qquad (29) $$

$$ \mathbf{C} = [0 \ 0 \ \cdots \ 0 \ 1], \qquad (30) $$

$$ 
\begin{aligned}
\mathbf{X}(k+1) = & \begin{bmatrix} \mathbf{A}_{x_a} & 0 & 0 & 0 \\ 0 & \mathbf{A}_{y_a} & 0 & 0 \\ 0 & 0 & \mathbf{A}_{x_b} & 0 \\ 0 & 0 & 0 & \mathbf{A}_{y_b} \end{bmatrix} \mathbf{X}(k) \\ 
& + \begin{bmatrix} \mathbf{B}_{x_a} & 0 & 0 & 0 \\ 0 & \mathbf{B}_{y_a} & 0 & 0 \\ 0 & 0 & \mathbf{B}_{x_b} & 0 \\ 0 & 0 & 0 & \mathbf{B}_{y_b} \end{bmatrix} \mathbf{u}(k),
\end{aligned} 
\qquad (31) $$

$$ \mathbf{Y}(k) = \begin{bmatrix} \mathbf{C}_{x_a} & 0 & 0 & 0 \\ 0 & \mathbf{C}_{y_a} & 0 & 0 \\ 0 & 0 & \mathbf{C}_{x_b} & 0 \\ 0 & 0 & 0 & \mathbf{C}_{y_b} \end{bmatrix} \mathbf{X}(k), \qquad (32) $$

where the linear part's MIMO system is expressed by rewriting the state equations (32) and (33); $\mathbf{X}(k)$ consists of

samples of the output $y$ before the $k$-th time step; $\mathbf{A}_{x_a}, \mathbf{A}_{x_b}, \mathbf{A}_{y_a}$, and $\mathbf{A}_{y_b}$ are composed of each motor model parameter matrix $\mathbf{A}$; similarly, other parameters of the MIMO system can also be expressed by equations (30) and (31).

### 3.2. Nonlinear Part.
The nonlinear part of the system mainly represents the relationship between the rotation of the puncture needle and the coordinates of the block. The relationship between them through equations (2)–(9) can be obtained as follows:

$$ (x_e, y_e, z_e) = \mathbf{f}(x_a, y_a, x_b, y_b), \qquad (33) $$

where $\mathbf{f}$ is their nonlinear mapping. Then, a first-order expansion of Taylor's formula is performed for $\mathbf{f}$ to obtain the following linear relationship:

$$ [x_e, y_e, z_e]^T = \mathbf{J}_f \begin{bmatrix} x_a - x_{ar} \\ y_a - y_{ar} \\ x_b - x_{br} \\ y_b - y_{br} \end{bmatrix} + \boldsymbol{\sigma}, \qquad (34) $$

$$ \mathbf{J}_f = \begin{bmatrix} \frac{\partial f_{x_t}}{\partial x_a} & \frac{\partial f_{x_t}}{\partial y_a} & \frac{\partial f_{x_t}}{\partial x_b} & \frac{\partial f_{x_t}}{\partial y_b} \\ \frac{\partial f_{y_t}}{\partial x_a} & \frac{\partial f_{y_t}}{\partial y_a} & \frac{\partial f_{y_t}}{\partial x_b} & \frac{\partial f_{y_t}}{\partial y_b} \\ \frac{\partial f_{z_t}}{\partial x_a} & \frac{\partial f_{z_t}}{\partial y_a} & \frac{\partial f_{z_t}}{\partial x_b} & \frac{\partial f_{z_t}}{\partial y_b} \end{bmatrix}, \qquad (35) $$

where $\mathbf{Y}_r(k) = [x_{ar}, y_{ar}, x_{br}, y_{br}]^T$ is the target reference values obtained by the inverse kinematics model, $\boldsymbol{\sigma}$ is an infinitesimal quantity, and $\mathbf{J}_f$ is the linearized Jacobian matrix. Finally, the complete input and output linear system model of the surgical robot can be obtained through equations (33) and (35) as follows:

$$ \mathbf{X}(k+1) = \mathbf{A}_s \mathbf{X}(k) + \mathbf{B}_s \mathbf{u}(k), \qquad (36) $$

$$ [x_e, y_e, z_e]^T = \mathbf{J}_f [\mathbf{C}_s \mathbf{X}(k) - \mathbf{Y}_r(k)], \qquad (37) $$

where

$$ \mathbf{A}_s = \begin{bmatrix} \mathbf{A}_{x_a} & 0 & 0 & 0 \\ 0 & \mathbf{A}_{y_a} & 0 & 0 \\ 0 & 0 & \mathbf{A}_{x_b} & 0 \\ 0 & 0 & 0 & \mathbf{A}_{y_b} \end{bmatrix}, \qquad (38) $$

$$ \mathbf{B}_s = \begin{bmatrix} \mathbf{B}_{x_a} & 0 & 0 & 0 \\ 0 & \mathbf{B}_{y_a} & 0 & 0 \\ 0 & 0 & \mathbf{B}_{x_b} & 0 \\ 0 & 0 & 0 & \mathbf{B}_{y_b} \end{bmatrix}, \qquad (39) $$
$$C_s = \begin{bmatrix} C_{x_a} & 0 & 0 & 0 \\ 0 & C_{y_a} & 0 & 0 \\ 0 & 0 & C_{x_b} & 0 \\ 0 & 0 & 0 & C_{y_b} \end{bmatrix}. \quad (40)$$

## 4. Model Predictive Control

Model predictive control (MPC) is now a well-known model-based control strategy, and it has gone through a period of development [41–43]. Taking the state-space model as an example, MPC uses the state of the model at the current time and the target state to estimate the control sequences at the next time under constraints. Compared with proportional-integral-derivative (PID) controller, MPC can predict the future behavior of the system. The optimal future control sequences are estimated by the known data. These sequences are evaluated at each sample time by the optimizer with constraints.

MPC consists of three components [44]: (1) a loss function, (2) constraints in the form of equality and inequality, and (3) initial conditions. The purpose of the loss function is to minimize the distance between the target state and the state at the current time [45], thereby minimizing the energy required by the control system. The equality constraints consist of the system’s state equations. The limits on the state variables are considered as the inequality constraints. If all constraints are linear, the problem forms a convex set. Meanwhile, if the control system is also convex, MPC can be solved as a convex optimization problem. Similarly, the MPC can be solved as a quadratic problem if the control system is a quadratic problem.

MPC uses the system’s model to predict the future control sequences [46]. As shown in Figure 5, the optimizer in MPC is to get the best control input. The linear time-invariant system will be considered in the discrete time domain as shown in equation (37). However, in equation (37), the state and input variables are added by two matrices to obtain the final output. To comprehensively consider the state and input variables, equations (37) and (38) are rewritten as follows:

$$\xi(k + 1) = \mathbf{A}_{\mathrm{m}}\xi(k) + \mathbf{B}_{\mathrm{m}}\Delta\mathbf{u}(k), \quad (41)$$

$$[x_e, y_e, z_e]^T + \mathbf{J}_f \mathbf{Y}_r(k) = \mathbf{C}_m \xi(k), \quad (42)$$

where

$$\xi(k) = \begin{bmatrix} \mathbf{X}(k) \\ \mathbf{u}(k-1) \end{bmatrix}, \quad (43)$$

$$\Delta\mathbf{u}(k) = \mathbf{u}(k) - \mathbf{u}(k-1), \quad (44)$$

$$\mathbf{A}_{\mathrm{m}} = \begin{bmatrix} \mathbf{A}_{\mathrm{s}} & \mathbf{B}_{\mathrm{s}} \\ 0 & I_{4\mathrm{nx}4} \end{bmatrix}, \quad (45)$$

$$\mathbf{B}_{\mathrm{m}} = \begin{bmatrix} \mathbf{B}_{\mathrm{s}} \\ I_{4\mathrm{nx}4} \end{bmatrix}, \quad (46)$$

$$\mathbf{C}_{\mathrm{m}} = [\mathbf{C}_{\mathrm{s}}, 0]. \quad (47)$$

Equations (42) and (43) are state-space equations that consider both state and input variables. We mainly discuss the optimization problem of equation (42), namely, the state space. This equation can make predictions about the future state, that is, the prediction step required for optimization can be obtained:

$$\mathbf{Y}_{\xi} = \mathbf{\Psi}_{\xi}(k) + \mathbf{\Theta}\Delta\mathbf{U}(k), \quad (48)$$

where

$$\mathbf{Y}_{\xi} = [\xi(k+1), \xi(k+2), \dots, \xi(k+N_p)]^T, \quad (49)$$

$$\mathbf{\Psi} = [\mathbf{A}_{\mathrm{m}}, \mathbf{A}_{\mathrm{m}}^2, \dots, \mathbf{A}_{\mathrm{m}}^{N_p}]^T, \quad (50)$$

$$\mathbf{\Theta} = \begin{bmatrix} \mathbf{B}_{\mathrm{m}} & 0 & 0 & \dots & 0 \\ \mathbf{A}_{\mathrm{m}}\mathbf{B}_{\mathrm{m}} & \mathbf{B}_{\mathrm{m}} & 0 & \dots & 0 \\ \mathbf{A}_{\mathrm{m}}^2\mathbf{B}_{\mathrm{m}} & \mathbf{A}_{\mathrm{m}}\mathbf{B}_{\mathrm{m}} & \mathbf{B}_{\mathrm{m}} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{A}_{\mathrm{m}}^{N_p-1}\mathbf{B}_{\mathrm{m}} & \mathbf{A}_{\mathrm{m}}^{N_p-2}\mathbf{B}_{\mathrm{m}} & \mathbf{A}_{\mathrm{m}}^{N_p-3}\mathbf{B}_{\mathrm{m}} & \dots & \mathbf{A}_{\mathrm{m}}^{N_p-N_c}\mathbf{B}_{\mathrm{m}} \end{bmatrix}, \quad (51)$$

where $N_p$ is the predicted step size of the system and $N_c$ is the control step size of the system. Next, the reference trajectory $\xi_r = [\mathbf{X}_r, \mathbf{u}_r]^T$ is used to construct the following loss function:

$$\mathbf{J}_{\xi} = [\mathbf{Y}_{\xi} - \mathbf{Y}_{\xi_r}]^T \mathbf{Q}_{\xi} [\mathbf{Y}_{\xi} - \mathbf{Y}_{\xi_r}], \quad (52)$$

where $\mathbf{Q}_{\xi}$ is the weight matrix. Considering that the control input cannot change too much within a period of time, the following loss function is constructed:

$$\mathbf{J}_{\mathrm{u}} = \Delta\mathbf{U}^T \mathbf{P}_{\mathrm{u}} \Delta\mathbf{U}, \quad (53)$$

where $\mathbf{P}_{\mathrm{u}}$ is the weight matrix. And the objective function to be optimized can be obtained by equations (53) and (54), which can be seen as a quadratic programming problem:

$$\begin{align*} \mathbf{J} &= \mathbf{J}_{\xi} + \mathbf{J}_{\mathrm{u}} \\ &= [\mathbf{\Psi}_{\xi} + \mathbf{\Theta}\Delta\mathbf{U}]^T \mathbf{Q}_{\xi} [\mathbf{\Psi}_{\xi} + \mathbf{\Theta}\Delta\mathbf{U}] - \mathbf{Y}_{\xi_r}^T \mathbf{Q}_{\xi} \mathbf{Y}_{\xi_r} + \Delta\mathbf{U}^T \mathbf{P}_{\mathrm{u}} \Delta\mathbf{U} \\ &= [\mathbf{\Psi}_{\xi}]^T \mathbf{Q}_{\xi} [\mathbf{\Psi}_{\xi}] - \mathbf{Y}_{\xi_r}^T \mathbf{Q}_{\xi} \mathbf{Y}_{\xi_r} + [\mathbf{\Theta}\Delta\mathbf{U}]^T \mathbf{Q}_{\xi} [\mathbf{\Theta}\Delta\mathbf{U}] \\ &\quad + 2[\mathbf{\Psi}_{\xi}]^T \mathbf{Q}_{\xi} [\mathbf{\Theta}\Delta\mathbf{U}] + \Delta\mathbf{U}^T \mathbf{P}_{\mathrm{u}} \Delta\mathbf{U} \\ &= \Delta\mathbf{U}^T [\mathbf{\Theta}\mathbf{Q}_{\xi}\mathbf{\Theta} + \mathbf{P}_{\mathrm{u}}] \Delta\mathbf{U} + 2[\mathbf{\Psi}_{\xi}]^T \mathbf{Q}_{\xi} [\mathbf{\Theta}\Delta\mathbf{U}]. \end{align*} \quad (54)$$

Finally, the constraints of the control variable $\Delta\mathbf{U}$ of the loss function to construct a standard form of quadratic programming will be set:
<figure><img src="image_6.png" /><figcaption>FIGURE 5: MPC controller internal rolling optimization process. The optimizer optimizes the current control input to obtain the optimal control input at the next step.</figcaption></figure>

$$ \begin{array}{ll} \min_{\Delta U} & J \\ \text{s.t.} & \text{lb} \le \Delta U \le \text{ub}, \end{array} \qquad (55) $$

where lb and ub are the lower and upper bounds of the variable $\Delta U$. Then, the best control sequences can be obtained by equation (45).

## 5. Simulation Setup and Results

This section mainly describes the simulated experimental settings and experimental results about this article. This article uses MATLAB to build a remote subsystem, which contains the MPC controller and parameter identification system. Under the function of the parameter identification system, the MPC controller obtains the robotic kinematics model realistically and outputs the optimal control sequence. In the first place, as shown in Figure 2, there is actually a TD when the MPC controller is applied to the system; thus, in the experiment, a delay system will be considered for comparison. We will, in addition, divide it into two experiments: (1) order and parameter identification of the system's linear part and (2) performance of MPC controller and PID controller in delay system.

### 5.1. Model Identification Experiment.
Above all, only the model identification of the linear part of the system is considered. And the order needs to be identified through offline data. The system's order can be obtained by equation (21). In this paper, the system's order $\hat{n} = 2$, the pitch of the screw $\rho = 0.02$, and the largest movement distance of the block $\eta = 2$. Next, the Kalman filter is designed to identify the system parameters. The state equation of the system is rewritten to identify the parameters of the system, and the result is represented by $\mathbf{H}$ in equation (28). The results of parameter identification are shown in Figure 6.

It can be seen from Figure 6 that the three parameters can stabilize at a certain value $a_1 = 0.0125$, $a_2 = -0.0257$, and $b = 1.818$ within the step size. However, it is not known whether a system composed of stable parameters can approximate the real system. Thus, as shown in Figure 7, the output of the actual system is compared with the output of the identified model. The figure shows the target curve, the output curve of the identification model, and the error curve

between them. The value of the error curve floats above and below zero.

### 5.2. MPC Performance.
First, the model of the linear part of the control system can be obtained by equation (28), and the parameters of the model are obtained in the previous experiment, which contains the effective length of the needle $h = 0.05$ m and the distance between blocks $d = 0.03$ m. And the linearized model of the nonlinear part of the control system can be obtained by equation (35). Then, the two models are combined by equation (37) to obtain a complete model of the control system. In addition, in order to comprehensively consider the control input and system state, the state-space equation of the control system is rewritten by equations (42) and (43). Finally, the prediction equation (49) required by the MPC controller is obtained by equation (42), which is used for rolling optimization in the MPC controller. The constraints of quadratic programming $lb = -0.002$ and $ub = 0.002$ are constraints on the difference values of control sequence. The optimizer of the MPC controller outputs the optimal control sequence through equation (55). In this paper, traditional PID controller is designed to compare with the MPC controller. In fact, as the controller takes time to run, there is a delay in the output of the controller, which is equivalent to the control delay system. Thus, we set different TDs in the control system. We set the weight matrix $Q_{\xi}$ and $P_u$ in the MPC controller as the identity matrix and the parameters of the PID controller $K_p = 13$, $K_i = 2$, and $K_d = 0$.

Figure 8 shows the Y-axis motion trajectory of block-a, which illustrates that with increase in the delay time, the response time of the system also increases and the desired trajectory can be tracked by the MPC controller. As the delay increases, the response speed of the system also becomes slower, and the output of the MPC controller will produce 2% overshoot. However, the PID controller does not perform well in the case of TD. The output of the PID controller diverges when TD = 250 ms, making the system unstable. When TD = 150 ms, the PID controller makes the system gradually stable, but the setting time ts = 4 s and the 10% overshoot exceed the constraints of the surgical robot.

Then, the output of the position of the tip of the puncture needle will be analyzed. As shown in Figure 9, in the open-loop state without controller, the output value of the system
<figure><img src="image_7.png" /><figcaption>FIGURE 6: The identification result of the system parameters, parameters $a_1$ and $a_2$ are the coefficients of the output term, and parameter $b$ is the coefficient of the input term.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7: Comparison of system parameter identification. The error is the difference between the target value and the estimated value.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8: The performance of the Y-axis motion trajectory of block A under different delays and different controllers.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 9: The performance of the Z-axis motion trajectory of the tip of puncture needle under different delays and different controllers.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10: The performance of the coordinate difference in the Y-axis direction of the block. The TD will cause the regular periodic error of the sine wave.</figcaption></figure>

has a large fluctuation around the desire trajectory. The closed-loop system formed by the MPC controller can track the desired trajectory with TD, which will produce 2% overshoot in the delay system. However, the closed-loop system formed by the PID controller is sensitive to TD. When there is a large TD in the system, it will produce a larger overshoot than the MPC controller and even make the system unstable.

Finally, in order to more clearly show the advantages of the MPC controller, a differential operation is performed on the state value of the sine wave shown in Figure 8 to obtain the distance from the target value shown in Figure 10. In the figure, the output of the MPC controller applied to the delay system produces phase shifts at different distances relative to the expected value. However,

when the PID controller acts on a nondelayed system, phase shift occurs. As the delay increases, even the output diverges. The above experimental results show that the effect of MPC controller is better than that of the PID controller. The error variance is shown in Table 1, the desired error variance is zero, and the closer to the ideal error variance, the more ideal the system's outputs are. The error variance of the MPC controller in the TD system is smaller than the error variance generated by the PID controller. And when TD = 0, the error variance generated by the MPC controller is even smaller than the error variance generated by the nondelay system with noise. From the result, we can get that for the parallel mechanism robot involved in this article, MPC has a superior control effect because MPC can control the MIMO system.
<table><caption>TABLE 1: Error variance (cm²) under different TD.</caption><thead><tr><th>Y1-MPC-250</th><th>Y1-MPC-150</th><th>Y1-MPC-50</th><th>Y1-MPC-0</th><th>Y1-pid-250</th><th>Y1-pid-150</th><th>Y1-pid-50</th><th>Y1-PID-0</th><th>Y1-reality-0</th></tr></thead><tbody><tr><td>0.6625</td><td>0.4048</td><td>0.1512</td><td>0.0728</td><td>1.6730</td><td>0.7067</td><td>0.6178</td><td>0.6860</td><td>0.1077</td></tr></tbody></table>

## 6. Conclusion and Future Work

We introduced a surgical robotic puncture needle positioning platform that can improve the efficiency of biopsy. Kinematic analysis of the platform illustrates that the puncture needle can work in the working space through multiple angles that can be obtained by decoupling the robotic kinematic. In this paper, two constraints are defined for decoupling the kinematic, which ensure the robotic largest workspace. The designed MPC controller makes the control flow of the positioning device form a closed loop. And the state equation of the MPC controller needed is divided into two parts: linear part and nonlinear part. In order to guarantee the simulation of robotic system closer to the reality, the Kalman filter system parameters identification approach is used to obtain the system's parameters in the linear part. In the nonlinear part, Taylor formula is used to linearize the system, and the linearization system only contains a Jacobian matrix. Next, the two parts of the system are combined into a complete robot system state equation. In simulation, the closed-loop control significantly improves the stability of robot needle positioning. And the performance of MPC controller is better than that of the PID controller. Finally, the system is open source in Github (https://github.com/tKsome/MPC-Puncture-Robot) to reduce the amount of development time for other researchers involved in the control of the medical robotic MIMO system. The system is flexible and may be useful for most image-based intervention procedures. In future work, we will test the performance of the proposed approach on the real-work robotic platform. Furthermore, a force control strategy will be proposed such that a force feedback can be provided to the surgeon through remote operation.

## Data Availability

The data included in this study are available upon proper request by contact with the corresponding author.

## Conflicts of Interest

The authors declare that they have no conflicts of interest.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China with grant no. 61803103 and Science and Technology Program of Guangdong with grant nos. 2019B010140002 and 2019B010142001.

## References

[1] P. Kulkarni, S. Sikander, P. Biswas, S. Frawley, and S.-E. Song, "Review of robotic needle guide systems for percutaneous

intervention," *Annals of Biomedical Engineering*, vol. 47, p. 7, 2019.

[2] B. Schell, K. Eichler, M. Mack et al., "Robotunterstützte punktion in einem hochfeld-kernspintomografen—erste klinische ergebnisse," *RöFo-Fortschritte auf dem Gebiet der Röntgenstrahlen und der Bildgebenden VERFAHREN*, vol. 184, no. 1, pp. 42–47, 2012.

[3] N. Hungr, I. Bricault, P. Cinquin, and C. Fouard, "Design and validation of a CT- and MRI-guided robot for percutaneous needle procedures," *IEEE Transactions on Robotics*, vol. 32, no. 4, pp. 973–987, 2016.

[4] N. Simaan, R. M. Yasin, and L. Wang, "Medical technologies and challenges of robot assisted minimally invasive intervention and diagnostics," 2018, https://arxiv.org/abs/1807.03731.

[5] Y. Chen, S. Xu, A. Squires et al., "MRI-guided robotically assisted focal laser ablation of the prostate using canine cadavers," *IEEE Transactions on Biomedical Engineering*, vol. 65, no. 7, pp. 1434–1442, 2017.

[6] S. Lim, K. Sharma, P. Li et al., "Robotically assisted long bone biopsy under MRI: cadaver study results," *International Journal of Computer-Assisted Radiology and Surgery*, vol. 14, no. 1, pp. 147–156, 2018.

[7] M. Giannakou, C. Yiallouras, G. Menikou, C. Ioannides, and C. Damianou, "MRI-guided frameless biopsy robotic system with the inclusion of unfocused ultrasound transducer for brain cancer ablation," *The International Journal of Medical Robotics and Computer Assisted Surgery: MRCAS*, vol. 15, p. e1951, 2019.

[8] B. Fallahi, C. Rossa, R. Sloboda, N. Usmani, and M. Tavakoli, "Partial estimation of needle tip orientation in generalized coordinates in ultrasound image-guided needle insertion," in *Proceedings of the 2016 IEEE International Conference on Advanced Intelligent Mechatronics (AIM)*, pp. 1604–1609, Banff, Canada, 2016.

[9] Y. Kurose, Y. M. Baek, Y. Kamei et al., "Preliminary study of needle tracking in a microsurgical robotic system for automated operations," in *Proceedings of the 2013 13th International Conference on Control, Automation and Systems (ICCAS 2013)*, pp. 627–630, Gwangju, Republic of Korea, 2013.

[10] A. Melzer, B. Gutmann, T. Remmele et al., "Innomotion for percutaneous image-guided interventions," *IEEE Engineering in Medicine and Biology Magazine*, vol. 27, no. 3, pp. 66–73, 2008.

[11] H. J. Won, N. Kim, G. Kim, J. B. Seo, and H. Kim, "Validation of a CT-guided intervention robot for biopsy and radio-frequency ablation: experimental study with an abdominal phantom," *Diagnostic and Interventional Radiology*, vol. 23, pp. 233–237, 2017.

[12] D. Navarro-Alarcon, S. Singh, T. Zhang et al., "Developing a compact robotic needle driver for MRI-guided breast biopsy in tight environments," *IEEE Robotics and Automation Letters*, vol. 2, no. 3, pp. 1648–1655, 2017.

[13] S. Park, J.-G. Kim, K.-W. Lim et al., "A magnetic resonance image-guided breast needle intervention robot system: overview and design considerations," *International Journal of Computer Assisted Radiology and Surgery*, vol. 12, pp. 1319–1331, 2017.
[14] T. Zhang, D. Navarro-Alarcon, K. W. Ng, M. K. Chow, Y. Liu, and H. L. Chung, "A novel palm-shape breast deformation robot for MRI-guided biopsy," in *Proceedings of the 2016 IEEE International Conference on Robotics and Biomimetics (ROBIO)*, pp. 527-532, Qingdao, China, December 2016.

[15] Y. Chen, A. Squires, R. Seifabadi et al., "Robotic system for MRI-guided focal laser ablation in the prostate," *IEEE/ASME Transactions on Mechatronics*, vol. 22, no. 1, pp. 107-114, 2017.

[16] Z. Cao, T. Liao, W. Song, Z. Chen, and C. Li, "Detecting the shuttlecock for a badminton robot: a yolo based approach," *Expert Systems with Applications*, vol. 164, Article ID 113833, 2020.

[17] W. Wang, A. N. Viswanathan, A. L. Damato et al., "Evaluation of an active magnetic resonance tracking system for interstitial brachytherapy," *Medical Physics*, vol. 42, no. 12, p. 7114, 2015.

[18] Y. Chen, W. Wang, E. J. Schmidt et al., "Design and fabrication of MR-tracked metallic stylet for gynecologic brachytherapy," *IEEE/ASME Transactions on Mechatronics*, vol. 21, no. 2, pp. 956-962, 2016.

[19] J. Ladoiye, D. Neculescu, and J. Sasiadek, "Force control of surgical robot with time delay using model predictive control," in *Proceedings of the 15th International Conference on Informatics in Control, Automation and Robotics*, Porto, Portugal, 2018.

[20] A. Squires, J. Oshinski, J. Lamanna, and Z. T. H. Tse, "Spiromet template: a platform for MRI-guided spinal cord injec-tions," *Journal of Medical Robotics Research*, vol. 1, no. 1, Article ID 1640006, 2016.

[21] L. Zhang, Z. Chen, W. Cui et al., "Wifi-based indoor robot positioning using deep fuzzy forests," *IEEE Internet of Things Journal*, 2020.

[22] W. Chen, T. Liao, Z. Li et al., "Using FTOC to track shuttlecock for the badminton robot," *Neurocomputing*, vol. 334, pp. 182-196, 2019.

[23] Y. Koseki, D. De Lorenzo, K. Chinzei, and A. Okamura, "Coaxial needle insertion assistant for epidural puncture," in *Proceedings of the 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems*, San Francisco, CA, USA, 2011.

[24] J. Kettenbach and G. Kronreif, "Robotic systems for percutaneous needle-guided interventions," *Minimally Invasive Therapy & Allied Technologies*, vol. 24, 2014.

[25] Y. Tipsuwan and M.-Y. Chow, "Control methodologies in networked control systems," *Control Engineering Practice*, vol. 11, no. 10, pp. 1099-1111, 2003.

[26] S. Zampieri, "Trends in networked control systems," *IFAC Proceedings Volumes (IFAC-PapersOnline)*, vol. 17, 2008.

[27] X.-M. Zhang, Q.-L. Han, X. Ge et al., "Networked control systems: a survey of trends and techniques," *IEEE/CAA Journal of Automatica Sinica*, vol. 7, no. 1, pp. 1-17, 2020.

[28] D. Uygun, "Low-power and cost-effective speed control card for remotely-controlled dc motors," *Electronics World*, vol. 121, pp. 16-19, 2015.

[29] T. Savi and M. Radonji, "One solution for remote control through internet application," in *Proceedings of the 2016 Information Technologies Conference*, Žabljak, Montenegro, 2016.

[30] H. Su, C. Yang, G. Ferrigno, and E. De Momi, "Improved human-robot collaborative control of redundant robot for teleoperated minimally invasive surgery," *IEEE Robotics and Automation Letters*, vol. 4, no. 2, pp. 1447-1453, 2019.

[31] H. Su, Y. Schmirander, S. Valderrama-Hincapie et al., "Asymmetric bimanual control of dual-arm serial

manipulator for robot-assisted minimally invasive surgeries," *Sensors and Materials*, vol. 32, p. 1223, 2020.

[32] C. Liu, W. H. Chen, and J. Andrews, "Model predictive control for autonomous helicopters with computational delay," in *Proceedings of the UKACC International Conference on Control*, Coventry, UK, 2010.

[33] P. Choomkasien, P. Chomphooyod, and D. Banjerdpongchai, "Design of model predictive control for industrial process with input time delay," in *Proceedings of the 2017 17th International Conference on Control, Automation and Systems (ICCAS)*, Jeju, Republic of Korea, October 2017.

[34] N. Dong, Y. Feng, X. Han, and A. Wu, "An improved model-free adaptive predictive control algorithm for nonlinear systems with large time delay," in *Proceedings of the 2018 IEEE 7th Data Driven Control and Learning Systems Conference (DDCLS)*, Enshi, China, May 2018.

[35] H. Norizuki and Y. Uchimura, "Design of model predictive control considering time delay in remote excavation work," in *Proceedings of the 2018 IEEE 15th International Workshop on Advanced Motion Control (AMC)*, Tokyo, Japan, March 2018.

[36] C. Cheng, S. Jiang, J. Liu, and J. Lou, "Design and kinematics analysis of a novel MR-compatible robot for needle insertion," "Design and kinematics analysis of a novel MR-compatible robot for needle insertion," in *Advances in Reconfigurable Mechanisms and Robots I*, J. S. Dai, M. Zoppi, and X. Kong, Eds., pp. 443-453, Springer, London, UK, 2012.

[37] S. Jiang, F. Sun, W. Feng, L. F. Hofman, and Y. Yu, "Analysis of a novel high-precision 5-degrees of freedom magnetic resonance imaging-compatible surgery robot for needle-insertion prostate brachytherapy," *Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science*, vol. 228, no. 5, pp. 865-876, 2014.

[38] E. G. Christoforou, I. Seimenis, E. Andreou, E. Eracleous, and N. V. Tsekos, "A novel, general-purpose, MR-compatible, manually actuated robotic manipulation system for minimally invasive interventions under direct MRI guidance," *The International Journal of Medical Robotics and Computer Assisted Surgery*, vol. 10, no. 1, pp. 22-34, 2014.

[39] S. Guo, Y. Wang, N. Xiao, Y. Zhao, Y. Zeng, and J. Wu, "A novel master-slave robotic system with close loop control for vascular interventional surgery," in *Proceedings of the 2018 IEEE International Conference on Mechatronics and Automation (ICMA)*, Tokyo, Japan, 2018.

[40] J. Guo, Y. Liu, Q. Qiu et al., "A novel robotic guidance system with eye gaze tracking control for needle based interventions," *IEEE Transactions on Cognitive and Developmental Systems*, p. 1, 2019.

[41] R. Kalman, "A new approach to linear filtering and prediction problems," *Journal of Basic Engineering*, vol. 82D, pp. 35-45, 1959.

[42] Y. Shang and T. Lu, "Optimal tracking control of mini unmanned helicopter based iterative approximating," in *Proceedings of the 2019 IEEE 3rd Information Technology, Networking, Electronic and Automation Control Conference (ITNEC)*, Chengdu, China, 2019.

[43] R. Tanaka and T. Koga, "An approach to linear active disturbance rejection controller design with a linear quadratic regulator for a non-minimum phase system," in *Proceedings of the 2019 Chinese Control Conference (CCC)*, Gangzhou, China, 2019.

[44] M. Dominici and R. Cortesão, "Model predictive control architectures with force feedback for robotic-assisted beating heart surgery," in *Proceedings of the 2014 IEEE International*
*Conference on Robotics and Automation (ICRA)*, Hong Kong, China, 2014.

[45] R. Abe, T. Satoh, N. Saito, J. Nagase, and N. Saga, "Disturbance observer-based model predictive control using time-varying constraints," in *Proceedings of the 2015 World Congress on Industrial Control Systems Security (WCICSS)*, London, UK, 2015.

[46] J. Jerez, E. Kerrigan, and G. Constantinides, "A condensed and sparse QP formulation for predictive control," in *Proceedings of the 2011 50th IEEE Conference on Decision and Control and European Control Conference*, Orlando, FL, USA, 2011.