Effects of parameters of particle swarm optimization algorithms on the path planning for flexible needles

Kinematic model

The kinematic model is the basis for path planning of flexible needles. Ideally, the trajectory of a flexible needle can be viewed as an arc with constant curvature, and the bending direction is determined by the orientation of the needle tip. An axial rotation of the needle can reorientate its tip. Thus, the movement of the flexible needle consists of a feed motion that performs the puncture and rotation motions changing the bending direction. That is, the control inputs for the flexible needle are insertion and rotation velocities (denoted by ${\mu }_1$ and ${\mu }_2$, respectively). Assuming that the motion of the needle tip determines that of the needle body and that rotations of the needle base are fully transmitted to the needle tip, the kinematic model can be described as follows

(1) $$\begin{array}{c}\dot{Q}=T\left(Q\right)\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right]\end{array}$$where Q represents the position and posture state of the needle tip in the world coordinate system and T(Q) is the transformation matrix between the control inputs and the needle tip states. These details can be found in Duindam et al. (2010), and Huo et al. (2018). As the radius of the needle trajectory is approximately independent of the insertion velocity (Webster, Memisevic & Okamura, 2005), and the motion of the needle base is perfectly transmitted to the needle tip (the above assumption), the insertion length of the needle (described as $l$) and rotation angle of the needle tip (represented as $\alpha$) can be used as the control inputs. Thus, the model can be discretized and rewritten as

(2) $$\begin{array}{c}{Q}_{i+1}=Q_i+T\left(Q_i\right)\left[\begin{array}{c}{l}_i\\ {\alpha }_i\end{array}\right]\end{array}$$where ${\alpha }_i$ and $l_i$ are the rotation angle and insertion length from the state $Q_i$ to $Q_{i+1}$, respectively. Based on the above analysis, a stop-and-rotate strategy (Wang et al., 2014) was chosen. For the ith arc, the operation involves axial rotation of the needle tip with angle ${\alpha }_i$ at the starting point $P_i$ and subsequent penetration over length $l_i$ to reach the end point $P_{i+1}$. Assuming that the initial puncture direction is perpendicular to the soft tissue, the intended path can be represented by the control sequence $\left[\left({\alpha }_1,l_1\right),\dots ,\left({\alpha }_i,l_i\right),\dots ,\left({\alpha }_n,l_n\right)\right]$.

To transform the control sequence into a planned path, inspired by the D-H method, arbitrary points on the planned path can be obtained by transforming the matrices (Wei & Jiang, 2021; Zhang et al., 2020). Figure 1 shows a schematic representation of this transformation. For the $ith$ arc, the coordinate $P_i{x}_i{y}_i{z}_i$ is attached to the needle tip following the right-hand rule. The $x_i$-axis points in the bending direction, and the $z_i$-axis is along the insertion direction perpendicular to the $x_i$-axis. The first operation to transform is to rotate the coordinate $P_i{x}_i{y}_i{z}_i$ around the $z_i$-axis with ${\alpha }_i$ angle. $Rot\left(z_i,{\alpha }_i\right)$ and $P_i{x}_i^{\prime }{y}_i^{\prime }{z}_i^{\prime }$ denotes this rotation operation and the obtained coordinate, respectively. Next, move $P_i{x}_i^{\prime }{y}_i^{\prime }{z}_i^{\prime }$ along the $x_i^{\prime }$-axis with length $a_i$ and along the $z_i^{\prime }$-axis with length $d_i$. The values of $a_i$ and $d_i$ can be easily calculated using geometry. After this translation operation denoted by $Trans\left(a_i,0,d_i\right)$, the coordinate becomes $P_{i+1}{x}_i^{″}{y}_i^{″}{z}_i^{″}$. Finally, rotate $P_{i+1}{x}_i^{″}{y}_i^{″}{z}_i^{″}$ around the $y_i^{″}$-axis by ${\theta }_i$. Here, ${\theta }_i$ represents the center angle of the $ith$ arc, which is obtained by dividing the radius r into arc length $l_i$. This rotation operation is represented as $Rot\left(y_i^{″},{\theta }_i\right)$. Thus, the transformation matrix from $P_i$ to $P_{i+1}$ can be expressed as

(3) $$\begin{array}{c}{T}_i^{i+1}\left({\alpha }_i,l_i\right)=Rot\left(z_i,{\alpha }_i\right)Trans\left(a_i,0,d_i\right)Rot\left(y_i^{″},{\theta }_i\right)\end{array}.$$

Assuming that the world coordinates coincide with those of the needle tip at the entry point, the position and pose of the needle tip after penetrating the $ith$ arc, represented as $T_1^{i+1}$, is described as follows:

(4) $$\begin{array}{c}{T}_1^{i+1}={\prod }_{j=1}^{j=i}{T}_j^{j+1}\left({\alpha }_j,l_j\right)\end{array}.$$

The states of the arbitrary points in the $ith$ arc can be calculated by replacing $l_j$ in Eq. (4) with the corresponding puncture length. Thus, a 3D kinematic model was created.

Path planning based on PSO

The basic idea of PSO is to obtain the optimal solution for a designed objective function by moving multiple particles (Cai et al., 2020). The key for path planning based on PSO is the design of particle. As the puncture length and rotation angle can be used for calculating the trajectory, the form of particle is $\left({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n;l_1,l_2,\dots ,l_n\right)$. These particles are updated according to Eqs. (5) to (6).

(5) $$\begin{array}{c}{v}_{i,g}=w{v}_{i,g-1}+c_1{r}_1\left(X_{ib}-X_{i,g-1}\right)+c_2{r}_2\left(X_{gb}-X_{ig-1}\right)\end{array}$$

(6) $$\begin{array}{c}{X}_{i,g}=X_{i,g-1}+v_{i,g}\end{array}$$where $v_{i,g}$ and $X_{i,g}$ represent the velocity and position of ith particle in gth generation, respectively. $w$, $c_1$ and $c_2$ are inertial weight coefficient and learning factors, respectively. They are the parameters needed to be optimized. $X_{ib}$ and $X_{gb}$ denote the best particles of ith particle and global, respectively. They are chosen according to the objective function. The objective function is determined considering the target deviation, path length, and obstacle avoidance. The path deviation $f_{pd}$ is expressed as follows:

(7) $$\begin{array}{c}{f}_{pd}=\left|P_{tip}-P_{target}\right|\end{array}$$where $P_{tip}$ and $P_{target}$ represent the endpoints of the given path and target, respectively. The path length $f_{pl}$ is expressed as follows:

(8) $$\begin{array}{c}{f}_{pl}=\sum _j^n{l}_j\end{array}.$$

During puncture, obstacle avoidance is pivotal because damage to sensitive tissues such as vessels and nerves can lead to severe complications. The path danger level $f_{pdl}$ is expressed as follows:

(9) $$f_{pdl}=\sum _i\{\begin{array}{cc}{\displaystyle \frac{DS}{d_i}}& \hfill if{d}_i>DS\\ +\mathrm{\infty }& \hfill otherwise\end{array}$$where $d_i$ is the minimum distance between the path and $ith$ barrier. DS is the safe distance between planned path and obstacles. The above factors are weighted differently ( ${\gamma }_1,{\gamma }_2,{\gamma }_3$), which are set to be 0.8, 0.1 and 0.1, respectively because we think that puncture accuracy is more important. Thus, the objective function is expressed as follows:

(10) $$\begin{array}{c}{f}_{of}={\gamma }_1{f}_{pd}+{\gamma }_2{f}_{pl}+{\gamma }_3{f}_{pdl}\end{array}.$$

Optimization by GA

The PSO parameters significantly influence the effects of path planning for flexible needles. As the GA is one of the most popular optimization algorithms with simple structure, it has been widely used in various applications. Thus, it is applied to obtain the best combination of parameters. GA evolutionarily determines the best individual using operators inspired by natural selection and biological processes. These operators include selection, crossover, and mutation. The form of GA chromosomes is designed as $\left(\omega ,c_1,c_2\right)$. Because the effects of the PSO parameters on path planning for flexible needles were studied through simulations based on the variable control method, a primary candidate combination of parameters was obtained. Thus, the primary combination was used as the best chromosome for the first generation, and the Lord strategy (Zhang et al., 2020) was applied instead of the roulette wheel selection, facilitating the convergence speed. When the Lord strategy is adopted, chromosomes are sorted according to the values of the fitness function, and the best one is selected as the Lord to perform a crossover with the chromosomes in even order. In our opinion, generally, the primary candidate combination obtained above performs much better than randomly generated chromosomes. Thus, selecting this solution rather than randomly generated one as the Lord helps to rapidly generate better chromosomes and facilitate the convergence speed. The crossover rate (CR) is empirically set to be 0.8. Note that the fitness is the average path deviation acquired by running the PSO 10 times with each specified combination of parameters. Subsequently, a multi-point mutation with mutation rate (MR) 0.1 was performed for each offspring. After the mutation, the parents and offspring are sorted and the number of population (NP) chromosomes are selected as the population of the next generation. Here, NP is 30. Note that NP here is the population size of GA. According to our experience, NP can be set to be 40 (Zhang et al., 2020). However, the computation time increases with the increasing of NP. Especially for optimizing the parameters of PSO, the PSO will be conducted about $10\times NP\times T$ times. T is the iteration times of GA. Thus, the computation time increases greatly with the increases of NP. To decrease the computation time, NP is set to be 30. This process was repeated 500 times (i.e., T is 500). These parameters are shown in Table 1. The optimization process is shown in Fig. 2.

Iterative PSO

As shown in ‘Optimization by GA’, even with the optimized combination of parameters, large deviations may occur between the path planned by PSO and the target. This may be caused by PSO falling into local minimum. To guarantee path planning accuracy, we propose the iterative PSO (I-PSO). The main idea of I-PSO is to iteratively running PSO until the path deviation within a specific accuracy (e.g., 0.001 mm). To reduce the computation time, parameters such as the number of particle swarms and iteration times are set according to the effects of PSO parameters. More details are shown in ‘Iterative PSO’.

Puncture platform

The established experimental puncture platform is shown in Fig. 3. A control system, two rotational stepping motors, and a digital camera are its primary components. The control system controls the stepping motors to perform the feeding and rotation movements. A digital camera with the resolution 2,592 × 1,944 was used to capture images during the puncture process to obtain the radius and trajectory of the needle insertion into the gelatin prostheses. The nitinol needle was customized using wire-electrode cutting. During the puncture verifying experiments, we first horizontally punctured the 18% gelatin using the flexible needle with 15° needle tip angle and 0.8 mm diameter at the speed of 15 mm/s for three times. Then, the puncture radius was calculated by image processing such as grey processing, denoising, sharpening and binarization. The puncture radius was used for path planning and the planned path would be punctured.

Data preprocessing is not involved in this experiment.

Effects of PSO parameters

First, the simulation environment was created (Scenario I, Fig. 4). The coordinates of the entry and target points were (0, 0, 0) and (100, 120, 120), respectively. The initial insertion direction was perpendicular to the x-y plane, and the radius was assumed to be 40 mm. The centers of the obstacles were (50, 60, 60) and (100, 90, 50) with a radius of 10 mm. A variable control method was then applied. As the learning factors, denoted by $c_1,c_2$, were set to 1 and 2 in Cai et al. (2020), Huo et al. (2018), Zhang et al. (2021), these values are adopted as the initial values of the PSO parameters to study the effects of inertial weight coefficient $\omega$. The population size was 50, and the number of iterations was 1,000. For each parameter combination, PSO was run five times and the average value will be chosen as the running results (shown in Fig. 5). Then simulations were done in MATLAB R2020a (The MathWorks, Natick, MA, USA) with AMD Ryzen 7 6800H and 16 G memory. The results show that $\omega$ has tiny effects on the computation time and small effects on the path length. However, when the value of $\omega$ changes from 0.1 to 1.5, the average path error varies from 0.9145 to 29.9272 mm, which significantly affects the puncture accuracy. For this condition, the best value for $\omega$ is 0.8.

Then, set the value of $\omega$ and $c_2$ to 0.8 and 2, respectively, and $c_1$ varies from 0.1 to 1.6 with an interval of 0.1. The results are shown in Fig. 6. From this, it can be seen that when $c_1$ increases, the average computation time is in the range of [2.0194, 2.1358], indicating that $c_1$ has little influence on the computation time. So as to effects of $c_1$ on the path length. However, for $c_1$, the average path error varied from 0.1585 to 10.2135 mm. The best $c_1$ value obtained under these conditions was 0.6.

Finally, under the condition $\omega =0.8,c_1=0.6$, the effects of $c_2$ is studied by simulations. The results are shown in Fig. 7. When $c_2$ changes in [0.1, 2.1], the average path error varies from 0.1585 to 12.5709 mm. The best value for $c_2$ in this scenario is 2. The above simulations demonstrate that $\omega$, $c_1$ and $c_2$ have little effect on the computation time of path planning for the flexible needle based on PSO and the path length, but significantly affect the path error, which is critical for needle insertion. To better show the effects of parameters, we take the effects of $c_1$ and $c_2$ when $\omega =0.8$ as an example. The results are shown in Fig. 8.

Optimization by GA

‘Effects of PSO parameters’ reveals a primary better combination of parameters of the form $c_1,c_2$ is [0.8, 0.6, 2]. This combination was used as the best-fit chromosome for the first-generation population of GA. After optimization, the best parameters were as $\left[0.6343,1.7289,1.8208\right]$. To illustrate the effectiveness of the optimized parameters, PSO was performed 100 times with these parameters. The results are shown in Fig. 9.

With the optimized $\omega ,c_1,c_2$, the path errors are in the range of [0, 11.0664]. More than half of the path errors are smaller than 0.1 mm (most are close to zero), but 35 path errors are still larger than 0.1 mm. This may be caused by the PSO falling into local minima. To overcome this problem, the effects of NP were studied by running PSO 100 times for each NP (Fig.10). Figure 10A show that with an increase of NP, the path error decreases rapidly at first and does not improve significantly when NP becomes larger than 100. Figure 10B shows the maximum path errors for different NP values. The NP value with the minimum average path error was 350. For NP = 350, the path error varied from 0 to 0.4886 mm, with an average of 0.0124 mm. This verifies the effectiveness of the parameter optimization. However, the computation time was approximately 12.6 s, which was unfavorable for intraoperative path replanning. Again, the simulation was based on MATLAB (The MathWorks, Natick, MA, USA) with a Ryzen 7 6800H and 16G memory.

Iterative PSO

To achieve better performance with less computation time, I-PSO was proposed based on the above results. The main idea is to execute the PSO iteratively until the desired accuracy is achieved. As described above, even with optimized parameters, PSO can plan a path with errors near zero even for a small population size. Because the computation time increases linearly with NP, a small NP helps reduce computation time. Thus, for I-PSO, NP was 50, T was 400, and the termination accuracy was 0.001 mm. One hundred times of simulations show that the path errors of the I-PSO are within [0, 0.00073] mm and 0.000046 mm on average, which significantly increases the planning accuracy compared to the optimized PSO with an average of 0.0124 mm. Moreover, the average computation time was 2.19 s, which was much shorter than that of the optimized PSO (12.6 s). The details are shown in Fig. 11.

Verifying experiments

To verify the performance of the PSO with the optimized parameters and the I-PSO, 100 times of simulation for each algorithm are executed and simulation results with Scenario I are compared with the existing simplified PSO (Cai et al., 2020) and adaptive intelligent PSO (AIPSO) (Tan et al., 2022a). For PSO, $\omega$, $c_1$ and $c_2$ are set to be 0.8, 1 and 2, respectively. For AIPSO, the parameters are the same with the reference, that is, ${\omega }_{max}=0.8$, ${\omega }_{min}=0.1$, $c_{max}=2$, $c_{min}=0.8$, $T=\mathrm{2,000}$ (Tan et al., 2022a). The results are shown in Table 2. From Table 2, we can summarize that the paths planned by optimized PSO are with least path length and smaller path error at the cost of increasing computation time. Paths planned by proposed I-PSO are with highest accuracy and the average computation time is acceptable. However, the path length is slightly increased. These results suggest that the proposed optimized PSO and I-PSO can significantly increase the path planning accuracy and compared to optimized PSO, I-PSO can greatly reduce the computation time (the average computation time of I-PSO is about 17% of that of optimized PSO).

To further verify the effectiveness of proposed algorithms, puncture experiments are carried out. Limited by the length of flexible needle and the depth of phantom carriage, we use the other scenario (Scenario II) for puncture. First, we establish the simulation environment and plan the path. Then, we execute the puncture experiment with the established puncture platform. The simulation environment of Cai et al. (2020) was used (Zhao et al., 2022). The target is located in (−5, 0, 120) and obstacle is located in (0, 0, 70) with a radius of 2 mm. The puncture radius was 407.24 mm, which was determined by a puncture experiment and image processing. The number of arcs is three. Because each time the path planned by PSO may be different, each algorithm is running only once for puncture experiment. According to literature, the path error of PSO for this scenario is 1.38 mm (Cai et al., 2020). To make the results in consistent with literature, $\omega$, $c_1$ and $c_2$ are set to be 1, 1 and 2, respectively. The simulation results show that, the optimized PSO can significantly improve the planning precision at the cost of computation time compared to PSO (Cai et al., 2020), and that I-PSO can achieve a sufficiently high precision with reduced time. The errors in the paths planned by both the optimized PSO and I-PSO are negligible. The path details are listed in Table 3.

The puncture experiments are shown in Fig. 12. During experiments, nitinol flexible needle with 0.8 mm diameter and 15 degrees needle tip is used to puncture the phantom with 18% gelatin at the speed of 15 mm/s. Due to the low permeability of gelatin, the target point and obstacle point in the experimental result diagram are not clear. Therefore, we mark the obstacle point and the target point, but do not modify the position of the target point and the obstacle point. At the same time, in order to better judge the position of the tip, the transparency of the target point mark is set to 40%. For the path planned by the optimized PSO, the coordinates of the needle tip for three times puncture are (−5.24, 0, 121.69), (−6.94, 0, 121.71) and (−7.97, 0, 120.11), which is 1.7070, 2.5861 and 2.9720 mm away from the ideal coordinates, respectively. Note that the coordinates are measured in the x-z plane; therefore, the y coordinate is set to 0, so that the error in the y-direction is neglected. For the path planned by I-PSO, the coordinate is (−5.96, 0, 121.09), (−6.21, 0, 120.6) and (−5.96, 0, 120.48), which is 1.4525, 1.3506 and 1.0733 mm off the ideal path, respectively. From Table 3 and puncture experiments, we can summarize that with the decreasing of path planning error, the puncture accuracy will be improved. The actual coordinates of the needle tip deviated from those of the simulations due to factors such as the simplification of the needle trajectory to an arc, out of sync between rotation angles of the needle base and tip and the accuracy of our puncture platform. The results also demonstrate that even with slight path planning error, the actual puncture error is much larger. The error of optimized PSO and I-PSO were within the acceptable range (4 mm according to Cai et al. (2020)).