Download PDF
Research Article  |  Open Access  |  21 Dec 2023

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Views: 285 |  Downloads: 60 |  Cited:   0
Intell Robot 2023;3(4):685-97.
10.20517/ir.2023.38 |  © The Author(s) 2023.
Author Information
Article Notes
Cite This Article


Effective management and monitoring of radioactive sources are crucial to ensuring nuclear safety, human health, and the ecological environment. A multi-robot collaborative radioactive source search algorithm based on particle swarm optimization particle filters is proposed. In this algorithm, each robot operates as a mobile observation platform using the latest observations to fuse into particle sampling. At the same time, the particle swarm optimization algorithm moves the particle set to a high-likelihood area to overcome particle degradation. In addition, each particle can learn from the search history of other particles to speed up the convergence of the algorithm. Lastly, the Dynamic Window Approach (DWA) for dynamic window obstacle avoidance is used to avoid obstacles in complex mountainous terrains to achieve efficient source search. Experimental results show that the search success rate of the proposed algorithm is as high as 95%, and its average search time is only 3.43 s.


Particle swarm optimization, particle filter, multi-robot, radioactive source search, DWA


In modern industrialized societies, nuclear technology is widely used in energy, medical treatment, scientific research, and other fields and has greatly contributed to human society. However, the International Atomic Energy Agency (IAEA) annual data report for 2022 revealed 146 incidents of illegal or unauthorized activities involving nuclear and radioactive materials [1]. Moreover, In March of 2023, a theft occurred in Salamanca, Guanajuato, Mexico, involving four sources of iridium-192, each with varying radioactivity levels of 35.64, 7.61, 1.14, and 0.11 Ci, along with the truck used for their transportation. As a result, an alert that covered seven central states was issued [2]. These activities pose serious consequences such as nuclear terrorism, radioactive contamination, and harm to individuals, presenting significant challenges to the global development of nuclear technology and nuclear safety.

Traditional source search methods for radioactive sources typically require significant human resources, which could be more efficient and pose considerable risks. In the early stages, robotic source detection involved using nuclear radiation detectors in a specific area, utilizing least-squares [3] and geometric methods [4] to estimate the location and intensity of the radiation source. For example, Howse et al. used recursive nonlinear least squares to estimate the location of the source, taking into account the variation in the source position and the deviation between the measurements and the model prediction [3]. Rao et al. used geometric differential triangulation to estimate the position and intensity of the source based on the measurements of three sensors [5]. In addition, the maximum likelihood estimation method is also widely used in source parameter estimation [6]. However, in the process of source parameter estimation, there is often the presence of noise and errors, which may affect the accuracy of the results. To address this issue, researchers have turned their attention to the application of Bayesian methods. This approach is based on Bayes' formula and combines prior knowledge and measurement data to calculate the posterior probability distribution of the radiation source, resulting in more accurate and reliable results [7, 8, 9, 10]. In addition, the Bayesian approach can utilize Markov chain Monte Carlo (MCMC) methods and particle filters (PF) to sample and compute the posterior probability distribution. Huo et al. employed a partially observable Markov decision process (POMDP) approach to enable autonomous robots to search for unknown radioactive sources [11]. Ling et al. introduced a collaborative search technique for unknown radioactive sources using multi-robot information fusion and an adaptive step-size free energy strategy [12]. In many practical scenarios, multi-robot systems face challenges in search tasks due to complex and limited feasible solutions. These challenges include effective resource utilization, maintenance of search capabilities, and avoidance of locally optimal solutions [13, 14]. Song et al. discussed the impact of different collaboration strategies and the number of robots on the required communication bandwidth and proposed a method to reduce communication bandwidth consumption by sharing the Gaussian mean and covariance matrix [15]. Effective multi-robot communication and data sharing are critical in these techniques since they allow the robots to exchange information on search tasks and achieve collaborative decision optimization.

Although many methods have been proposed to solve the problem of radioactive source localization in unknown areas, it remains a challenging task. The main difficulties come from: (1) The traditional single-robot algorithm has low resource utilization and small search coverage; (2) Due to the performance limitations of the detector and the influence of environmental factors, noise and errors may occur in the source parameter estimation process. Therefore, this paper proposes a multi-robot collaborative search algorithm based on particle swarm optimization (PSO)-PF [Figure 1]. The main contributions of this paper are reflected in the following points: (1) An improved PSO-PF algorithm based on multi-robots is proposed. This algorithm effectively finds the minimum difference between the latest observation and prediction results of each robot in the particle collection by minimizing the fitness function, thereby improving the accuracy of the PF; (2) The algorithm is able to find lost radioactive sources while avoiding obstacles in complex mountainous terrains, extending its applicability and practicality. The remainder of this paper is organized as follows. Section 2 outlines the radioactive source positioning task and provides a comprehensive explanation of the PF algorithm, including its principles and mathematical derivations. Section 3 delves into the specifics of the PSO-PF algorithm and the multi-robot path planning strategy proposed in this study. Section 4 presents simulation experiments and conducts an in-depth parameter analysis. Finally, Section 5 summarizes the key findings of this research and explores potential avenues for future research.

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Figure 1. The outline of the proposed algorithm architecture.


2.1. Problem description

Assume that the missing radioactive source $$ \theta=\left[x_{s}, y_{s}, I_{s}\right]^{T} $$ is located on the ground level in a mountainous area, where $$ \left(x_{s}, y_{s}\right) $$ is the position coordinate of the radioactive source, and $$ I_{s} $$ represents the activity of the source. In this study, a team of mobile robots was deployed to search for a lost radioactive source $$ \theta $$ in a complex mountainous environment. Each robot has the ability to move autonomously and perform local search operations while seamlessly cooperating with other robots to achieve a comprehensive and systematic global search.

Nuclear radiation decay is a random process that occurs in radionuclides. This randomness makes the decaying particles uncertain. According to previous studies, radiation counts from nuclear decay obey Poisson statistics [16]. Therefore, the probability that a radiation detector registers $$ z \in Z^{+} $$counts per second (CPS) from the source that emits on average $$ m $$ CPS is:

$$ \begin{equation} \mathcal{P}(z ; \lambda(\theta))=\frac{\lambda^{z}}{z !} e^{-\lambda} \end{equation} $$

where $$ \lambda=m \varepsilon $$ denotes the actual average count of particles detected by the detector per second and the parameter of the Poisson distribution; $$ \varepsilon $$ is the intrinsic detection efficiency of the detector, which is the ratio of the number of particles detected by the detector to the total number of particles incident on the detector at the same time. $$ \lambda=m \varepsilon $$ represents the true average count measured by the detector per second; $$ \varepsilon $$ is the inherent detection efficiency of the detector, which is the ratio of the number of particles detected by the detector to the total number of particles incident on the detector at the same time. Given the source parameter vector $$ {\theta_k} $$, the likelihood function of the radiation count measurement $$ z_{k} $$ at the recorded measurement location $$ \left(x_{k}, y_{k}\right), k=1, 2, \ldots $$ is given by:

$$ \begin{equation} p\left(z_{k} \mid {\theta_k}\right)=\mathcal{P}\left(z_{k} ; \lambda( {\theta_k})\right) \end{equation} $$

where $$ \mathcal{P} $$ is the Poisson distribution defined in Equation (1), and $$ \lambda( {\theta_k}) $$ is the mean radiation count recorded by the detector. When the distance between the radiation detector and the radiation source is $$ r $$, the calculation of $$ \lambda( {\theta_k}) $$ is performed as follows:

$$ \begin{equation} \lambda( {\theta_k})=\frac{I_{s} \times e^{\sum_{i}-u_{i} R_{i}}}{4 \pi r^{2}}+I_{b} \end{equation} $$

where $$ R_{i} $$ and $$ u_{i} $$ are the distance and radiation attenuation coefficient of the gamma rays through the ith medium, respectively; $$ I_{b} $$ is the background radiation count, which is the average count per second when $$ I_{s}=0 $$ or when $$ r $$ is infinite. The radiation detector measures the intensity of the radiation source in CPS and calculates the count value. Simultaneously, the energy response constant $$ \rho\left(C P S \cdot \mu S v^{-1} \cdot h\right) $$ is considered [17]. The count value $$ z $$ is expressed as the product of the mean radiation count $$ \lambda( {\theta_k}) $$ and the energy response constant $$ z=\lambda( {\theta_k}) \times \rho $$.

The activity of a radioactive source can be greatly affected by factors such as detector performance limitations, environmental absorption, multipath attenuation, and reflections. These factors often bring errors to detection data, thus affecting the accuracy of radioactive source positioning and parameter estimation. Figure 2 depicts the radiation profile attenuated by simulated radiation data according to Equation(3). It can be clearly observed from Figure 2 that due to the influence of the object's radiation attenuation coefficient, the radiation of the radioactive source model presents a complex distribution.

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Figure 2. Heatmap of the radiation distribution in the simulated environment, where the color bars indicate the signal strength in CPS. In the map, the source is at the upper right with an activity of about $$ 2.8 \times 10^5 $$ CPS. CPS: Counts per second.

2.2. Particle filter

Within the framework of the PF algorithm, we use an ensemble of particles and their associated particle weights to estimate the posterior probability density of the system state. Each particle represents a potential system state, and the accuracy of state estimation increases with the number of particles. At the same time, each particle weight represents the likelihood probability associated with the system state of each particle. Compared with traditional Kalman filtering methods, particle filtering has proven to be more suitable for dealing with nonlinear and non-Gaussian systems [18].

Suppose we have $$ N $$ mobile robots deployed in a radiation field, each maintaining $$ M $$ random particles. The estimated source parameter vector of the $$ m $$ th particle after obtaining the $$ k $$ th observation is expressed as $$ \theta_{k}^{m}=\left[x_{k}^{m}, y_{k}^{m}, I_{k}^{m}\right]^{T}, m \in M $$. At the $$ k $$ th moment, the $$ i $$ th robot is located at the location $$ \left(x_{k}^{i}, y_{k}^{i}\right) $$, and the obtained observation vector is expressed as $$ g_{k}^{i}= $$$$ \left\lceil x_{k}^{i}, y_{k}^{i}, z_{k}^{i}\right\rceil^{T}, i \in N $$. If the initial confidence is given as $$ p\left(\theta_{0} \mid g_{0}\right)=p\left(\theta_{0}\right) $$, then according to the Markov hypothesis, the prior probability density can be expressed as:

$$ \begin{equation} p\left(\theta_{k} \mid g_{1: k-1}\right)=\int p\left(\theta_{k} \mid \theta_{k-1}\right) p\left(\theta_{k-1} \mid g_{1: k-1}\right) d \theta_{k-1} \end{equation} $$

where $$ P\left(\theta_{k} \mid \theta_{k-1}\right) $$ represents the state transition probability, $$ p\left(\theta_{k-1} \mid g_{1: k-1}\right) $$ is the posterior probability density corresponding to the previous step. After obtaining the observation vector $$ g_{k} $$, the prior value is updated through Bayesian theory, and the posterior probability density is obtained as:

$$ \begin{equation} p\left(\theta_{k} \mid g_{1: k}\right)=\frac{p\left(g_{k} \mid \theta_{k}\right) p\left(\theta_{k} \mid g_{1: k-1}\right)}{p\left(g_{k} \mid g_{1: k-1}\right)} \end{equation} $$

where $$ p\left(g_{k} \mid g_{1: k-1}\right)=\int P\left(g_{k} \mid \theta_{k}\right) P\left(\theta_{k} \mid g_{1: k-1}\right) d \theta_{k} $$ is a normalization constant, which depends on the likelihood function $$ p\left(g_{k} \mid \theta_{k}\right) $$ and the statistical properties of the measurement noise. In the PF algorithm, we use the Monte Carlo (MC) method to approximate the posterior probability density of the target system state by the weighted sum of a set of particles [9].

$$ \begin{equation} p\left(\theta_{k} \mid g_{1: k}\right) \approx \sum\limits_{m=1}^{M} w_{k}^{m} \delta\left(\theta_{k}-\theta_{k}^{m}\right) \end{equation} $$

where $$ M $$ is the number of particles, $$ w_{k}^{m} $$ is the weight of the $$ m $$ th particle, and $$ \delta(\cdot) $$ is the Dirac function. Furthermore, the importance sampling (IS) technique [19] and sequential IS (SIS) [18] technique can be used to approximate the weights of these random particles.

$$ \begin{equation} w_{k}^{m} = w_{k-1}^{m} \cdot p\left(g_{k} \mid \theta_{k}^{m}\right) \end{equation} $$

where $$ p(\cdot) $$ is the target distribution, representing the posterior probability density we hope to obtain. In the context of radiation source localization, the fusion of measurements from multiple sensors is crucial for achieving accurate estimations of source parameters. Consider a scenario where an array of sensors is deployed to monitor radiation levels using a total number of sensors. To effectively use these measurements, we formulate the estimation of the radiation source parameters as a probabilistic exhibit minimal cross-correlation in their measurement errors. The unnormalized weight for the $$ i $$ th potential parameter value at time $$ k $$. is expressed as:

$$ \begin{equation} {w^{\prime}}_{(k)}^{(i)}=\prod\limits_{j=1}^{N} p\left(g_{(k, j)} \mid \theta_{(k-1)}^{(i)}\right) \cdot w_{(k-1)}^{(i)} \end{equation} $$

where $$ {w^{\prime}}_{(k)}^{(i)} $$ denotes the unnormalized weight associated with the approximated posterior probability density function. Consequently, the normalized weight for the $$ i $$ th potential parameter value is given by:

$$ \begin{equation} w_{(k)}^{(i)}=\frac{w_{(k)}^{\prime(i)}}{\sum_{j=1}^{N} w_{(k)}^{(j)}} \end{equation} $$

The PF is a non-parametric filtering algorithm extensively employed for state estimation tasks. Nevertheless, a major drawback of the PF is the issue of particle degeneracy (a few particles end up having non-zero weights). This degeneracy occurs when the likelihood function becomes highly localized, indicating precise observation information. Consequently, the overlap between the likelihood probability distribution and the prior probability distribution diminishes significantly, resulting in an inadequate number of particles to represent the posterior distribution effectively. Under such circumstances, the estimation outcome of the PF may exhibit considerable errors or even yield incorrect estimations. To address the issue of particle poverty, it becomes imperative to employ a larger number of particles to approximate the posterior distribution more accurately. However, this approach comes at the cost of increased computational complexity for the algorithm. Moreover, the effectiveness of the PF is contingent upon the accuracy and adaptability of the proposed distribution. When a substantial disparity exists between the proposed distribution and the target distribution, the performance of the PF is adversely affected. For instance, if the target distribution exhibits multiple peaks and the proposed distribution fails to accurately capture these peaks, the PF performance will be compromised. Furthermore, the PF also has the problem of high computational complexity. In practical applications, a large number of particles are often required to approximate the posterior distribution adequately. However, as the number of particles increases, so does the computational burden of the algorithm. This limitation restricts the real-time applicability of PFs in certain scenarios.


This section introduces a novel method that employs an improved PSO algorithm based on multi-robots to address the limitations of traditional PF algorithms, such as particle poverty and reliance on proposal distributions. This algorithm makes full use of the global search capability of PSO and the state estimation capability of PF algorithms to achieve highly accurate state estimation results within a relatively short period of time. Furthermore, the algorithm demonstrates remarkable robustness and adaptability, making it suitable for various system state estimation problems.

3.1. Improved particle swarm optimization algorithm based on multi-robot

PSO is an advanced global optimization algorithm that draws inspiration from the collective behavior of organisms such as birds or fish. It was originally introduced by Kennedy and Eberhart in 1995 [20]. This algorithm efficiently explores multi-dimensional search spaces by simulating the movement and interaction of a swarm of particles. In PSO, every potential solution within the search space assumes the form of a particle, possessing both a position and a velocity. As the particles diffuse, each particle adjusts its position and velocity based on its current state while retaining knowledge of its own historical best position and the overall best position found so far. Through continuous iterative refinement, the ultimate global optimal solution gradually reveals itself. Where the velocity and position of each particle are updated using the following formula:

$$ \begin{equation} \left\{\begin{array}{l} v_{t+1}^{i}=w v_{t}^{i}+c_{1} r_{1}\left(\text { pbest }^{i}-x_{t}^{i}\right)+c_{2} r_{2}\left(\text { gbest }-x_{t}^{i}\right) \\ x_{t+1}^{i}=x_{t}^{i}+v_{t+1}^{i} \end{array}\right. \end{equation} $$

where $$ x_{t}^{i} $$ signifies the position of the ith particle at time t; $$ v_{t}^{i} $$ symbolizes the velocity of the ith particle at time t; $$ \text{pbest}^{i} $$ reflects the best historical position reached by the ith particle, and gbest represents the best global position discovered so far. $$ w $$ is the inertia weight, $$ c_{1} $$ and $$ c_{2} $$ are learning factors, and $$ r_{1} $$ and $$ r_{2} $$ are random numbers between 0 and 1.

The algorithm presented in this paper integrates the PSO algorithm into the PF algorithm, leveraging it as the generator of the proposed distribution. The algorithm uses an adaptive weighted PSO algorithm, which dynamically adjusts inertia weights and acceleration coefficients to strike a balance between global and local search capabilities. The adaptive weight of the algorithm can be calculated as follows:

$$ \begin{equation} w(t)=w_{\min }-\frac{\left(w_{\max }-w_{\min }\right) t}{T} \end{equation} $$

where $$ w_{\max } $$ is the maximum value of the inertia weight, $$ w_{\min } $$ is the minimum value, and $$ T $$ denotes the maximum number of iterations. Furthermore, an adaptive fitness function is introduced to fuse the latest observations and recommendation distributions from multiple robots, defined as follows:

$$ \begin{equation} \varGamma =\sum\limits_{i=1}^{n} a b s\left(y_{\text {new }}^{i}-y_{\text {pre }}^{i}\right) \end{equation} $$

where $$ y_{\text {new }}^{i} $$ and $$ y_{\text {pre }}^{i} $$ denote the most recent observed and predicted values of the j-th robot, respectively; $$ \varGamma $$ represents the absolute deviation value between the observed value of the robot and the predicted value of the particle set at time $$ k $$. This characteristic is evident in the likelihood function $$ p\left(z_{k} \mid {\theta_k}\right) $$ for the radiation measurement value $$ y_{\text {new }}^{i} $$, wherein the weight increases as the measurement value approaches the peak and decreases as it deviates away from the peak. This paper presents an improved PSO algorithm based on multi-robot systems as a comparison to the conventional PSO method [2123]. The proposed algorithm incorporates the most recent observation and prediction outcomes from every robot, allowing the particles of all robots to relocate towards regions with higher posterior probability density distributions. This adjustment aims to overcome the issue of particle degradation, which occurs when the weight distribution of particles becomes highly imbalanced. Algorithm 1 outlines the process of particle optimization using the multi-robot improved PSO algorithm.

Algorithm 1 Improved Particle Swarm Optimization (IPSO)
    Input: Fitness function obj_fun, Lower bound lb, Upper bound ub, Population size pop_size, Current detector xk, Current particles cur_particles, Velocities vel, Personal best positions pbest, Global best position gbest, Velocity limit v_limit
    Output: Global best position gbest, Personal best positions pbest, Particle positions pos, Velocities vel, Iteration data iter_data
1: $$ W_{\text{max}}, c_{1}, c_{2}, \text{max_iter}, W_{\text{min}} \gets $$ Initialize variables
2: pos$$ \gets $$ Extract positions from cur_particles
3: for$$ i \gets 1 $$ to max_iterdo
4:        $$ W_{\text{pso}} \gets $$ Update inertia weight using Equation (11)
5:        for$$ n \gets 1 $$ to pop_sizedo
6:                v_id$$ \gets $$ Update velocity using Equation (10)
7:                Clip v_id to the range [-v_limit, v_limit]
8:                x_id$$ \gets$$x_id + v_id
9:                if (x_id(1) > ub and v_id(1) > lb) or (x_id(1) < lb and v_id(1) < 0) then
10:                    Reflect on the boundary for the $$ x $$-direction
11:                end if
12:                if (x_id(2) > ub and v_id(2) > lb) or (x_id(2) < lb and v_id(2) < 0) then
13:                    Reflect on the boundary for the $$ y $$-direction
14:                end if
15:                Update pos and vel for particle $$ n $$
16:                Compute fitness values $$ w_{\text{particles}} $$, $$ w_{\text{pbest}} $$, and $$ w_{\text{gbest}} $$
17:                if$$ w_{\text{particles}} < w_{\text{pbest}} $$then
18:                    Update personal best for particle $$ n $$
19:                end if
20:                if$$ w_{\text{particles}} < w_{\text{gbest}} $$then
21:                    Update global best
22:                end if
23:        end for
24:        if$$ w_{\text{pbest}} \times 0.9 \leq w_{\text{gbest}} $$then
25:            Break the loop if convergence criteria are met
26:        end if
27: end for
28: iter_data$$ \gets $$ Array containing iteration and global best fitness
29: returngbest, pbest, pos, vel, iter_data

3.2. Dynamic Window Approach dynamic window avoidance strategy

The Dynamic Window Approach (DWA) algorithm is a widely used method for obstacle avoidance [24]. Assuming the current state of the robot is denoted as $$ x_{t}=[x, y, \theta, v, \omega]^{T} $$, where $$ \mathrm{x} $$ and $$ \mathrm{y} $$ represent the coordinates of the robot, $$ \theta $$ represents the orientation, and $$ \mathrm{v} $$ and $$ \omega $$ represent the linear and angular velocities, respectively. The DWA algorithm generates a set of speed commands, which can be expressed as follows:

$$ \begin{equation} \left\{\begin{array}{l} v_{m}=\left\{(v, w) \mid v \in\left[v_{\text {min }}, v_{\text {max }}\right] \cap w \in\left[w_{\text {min }}, w_{\text {max }}\right]\right\} \\ v_{d}=\left\{(v, w) \mid v \in\left[v_{a}-\dot{v} * \Delta t, v_{a}+\dot{v} * \Delta t\right] \cap w \in\left[w_{a}-\dot{w} * \Delta t, w_{a}+\dot{w} * \Delta t\right]\right\} \\ v_{a}=\left\{(v, w) \mid v \in\left[v_{\text {min }}, \sqrt{2 * \operatorname{dist}(v, w) * \dot{v}}\right] \cap w \in\left[w_{\min }, \sqrt{2 * \operatorname{dist}(v, w) * \dot{w}}\right]\right\} \end{array}\right. \end{equation} $$

where $$ v_{\min } $$ and $$ v_{\max } $$ represent the lower limit and upper limit of the robot's linear speed, respectively; $$ w_{\min } $$ and $$ w_{\max } $$ present the lower limit and upper limit of the robot's angular speed, respectively; $$ v_{a} $$ and $$ w_{a} $$ represent the current linear speed and angular velocity of the robot, respectively; $$ \operatorname{dist}(v, w) $$ represents the shortest distance between the simulated trajectory and the obstacle at the current speed; $$ \Delta t $$ represents the time interval required for the robot to execute the speed command.

The DWA algorithm generates multiple sets of feasible speed command combinations, represented as $$ u=[v, w]^{T} $$. These combinations are determined by the intersection of velocity boundary limits $$ v_{m} $$, acceleration limits $$ v_{d} $$, and environmental obstacle limits $$ v_{a} $$. Figure 3 visually illustrates the predicted trajectories derived from these combinations of velocity commands using the kinematic model of the mobile robot. In addition, weights are assigned to these trajectories through the evaluation function $$ g(u) $$ to select the driving speed corresponding to the optimal trajectory.

$$ \begin{equation} g(u)=w_{1} \cdot \operatorname{dist}(u)+w_{2} \cdot \operatorname{heading}(u)+w_{3} \cdot \operatorname{speed}(u) \end{equation} $$

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Figure 3. Mountain multi-robot trajectory prediction model. The mobile robot's historical search trajectories are depicted in different colors: blue, orange, yellow, and purple dots. The solid green line represents the predicted trajectory, while the green triangles along the trajectory represent the expected robot position.

where $$ \operatorname{dist}(u) $$, heading $$ (u) $$, and speed $$ (u) $$ represent the distance, angle, and speed, respectively, between the robot's next position $$ x_{k+1} $$ and the target position under decision $$ u $$. The weight coefficients, $$ w_{1}, w_{2} $$, and $$ w_{3} $$, determine each factor's relative importance in the speed command evaluation function. When a trajectory intersects with an obstacle at any point in the future, that trajectory is disregarded. Finally, the optimal decision $$ u^{*} $$ is decided to minimize $$ g(u) $$ as:

$$ \begin{equation} u^{*}=\arg \min\limits_{u \in U} g(u) \end{equation} $$

where the feasible decision set $$ U $$ refers to the set of permissible speed instructions, as defined by Equation (13). Algorithm 2 outlines the process of computing the reward function for the nth robot at time $$ k $$.

Algorithm 2 Calculating the reward function for the nth robot at time $$ k $$.
    Input: Robot state matrix x, Estimated location matrix $$ goal $$
    Output: Robot decision matrix $$ u_{k} $$, Reward value matrix $$ G(u_{k}) $$
1: $$ mode $$$$ \gets $$ set kinematic limit parameter matrix
2: $$ \alpha , \gamma \longleftarrow $$ set the percentage constants for the evaluation
3: $$ \left\{mode, x_{k}\right\} \rightarrow \left\{v_m, v_d\right\} $$ feasible decision sets using Equation (13)
4: for all $$ u_{k}\in\{v_m, v_d\} $$do
5:        $$ \{u_{k}, x_{k}\} \rightarrow \{x_{k+1}, traj\} $$ Generate Trajectory
6:        $$ \{dist(u), heading(u), speed(u)\} \longleftarrow \{x_{k+1}, traj, goal\} $$ using Equation (15)
7:        for i = 1, ..., length(traj) do$$ \leftarrow $$ Iterate through all traj
8:            if$$ traj(i) $$ intersects with an obstacle then
9:                isIntersecting = true
10:            end if
11:        end for
12:        if$$ isIntersecting = true $$then
13:            continue;
14:        end if
15:end for
16: $$ \{w_1, w_2, w_3\} \rightarrow G(u_{k}) $$ using Eq.(18)
17: return$$ \{u, G\} $$ returns all decision Matrix u and its values G


4.1. An illustrative run

The simulation results of multi-robot collaborative radioactive source search using PSO-PF algorithms are shown in Figure 4. The illustrative run uses the following parameters for the simulation:

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Figure 4. An illustrative run of the source searching process using a multi-robot system. (A-C) Trajectories of the four robots and the estimated source locations at steps $$ k=1 $$, $$ k=5 $$, and $$ k=18 $$, respectively; (D) the PSO iteration and global fitness; (E) the observables of agents during the search. PSO: Particle swarm optimization.

● Multiple robots search for lost radioactive sources in a $$ 150 \mathrm{\; m} \times 150 \mathrm{\; m} $$ mountainous environment.

● The radioactive source is located at coordinates $$ (x_{s}=135 \mathrm{\; m}, y_{s}=90 \mathrm{\; m}) $$, and its activity is set to $$ I_{s}=2.94 \times 10^{6} \mathrm{\; CPS} $$.

● The intrinsic detection efficiency of the radiation detector $$ \varepsilon=60\% $$, the mountain radiation attenuation coefficient $$ u=0.05 $$, the energy response constant $$ \rho = 200 \mathrm{\; CPS/\mu Sv/h} $$, and the background radiation count is $$ I_{b}=20 \mathrm{\; s}^{-1} $$.

● The robot's movement is constrained by various physical parameters, including maximum linear velocity $$ v_{\max }=3 \mathrm{\; m/s} $$, maximum angular velocity $$ w_{\max }=90^{\circ}/\mathrm{s} $$, linear acceleration $$ \dot{v}=1 \mathrm{\; m/s}^{2} $$, angular acceleration $$ \dot{w}=30^{\circ}/\mathrm{s}^{2} $$, linear velocity resolution $$ 0.5 \mathrm{\; m/s} $$, and angular velocity resolution $$ 15^{\circ}/\mathrm{s} $$.

Figure 4A-C visually shows the source search area, the search trajectories of the robot at times $$ k $$ = 1, 5 and 18 using this algorithm, and the contour plots of the radiation distribution from the source. It is worth noting that the robot's trajectory, position, observation values, and particles are represented by different colors. Furthermore, the observations are described in the form of a staircase plot, and the global fitness and number of iterations are represented as a stacked plot of two variables with a common $$ \mathrm{x} $$-axis. During the search process, the fitness function of the PSO algorithm is calculated based on the latest observation values obtained by each robot, which makes the particle set move to a high-likelihood area and overcomes particle degradation. At step $$ k $$ = 1, the measurements obtained by the detector are small, and the positions of the particles are evenly distributed on the map. As the particle swarm algorithm iterates to find the optimal particle set, the particles gather around the estimated source location at step $$ k $$ = 5, indicating that the distance between the estimated source location and the true source location is decreasing. Finally, all robots complete the source search at step $$ k $$ = 18.

Figure 4D and E shows that as the number of search steps increases, the robot's observations become increasingly larger, and the number of iterations of the PSO algorithm gradually decreases. It should be emphasized that when the robots are far away from the source location, the observations obtained by the robots from each other tend to be relatively small. Therefore, as the difference between the observations obtained by the robot and its predicted observations becomes smaller, the global optimal particle fitness value also becomes smaller. On the contrary, when the robot approaches the source position and obtains relatively large observed values, the global optimal particle obtains higher fitness values that increase as the gap between the observed and predicted values increases.

4.2. Parametric analysis

In this experiment, we aim to evaluate the performance comparison of the proposed PSO-PF algorithm and the Extended KalmanParticle Filter algorithm (EPF) in the multi-robot radioactive source search problem. The termination conditions for the source search process are defined as follows: The process concludes when the Euclidean distance between all robots and the estimated source position is less than 5 m and the variance of the samples in the PF is less than 1. Once the search process ends, it is deemed successful if the estimated source locations have converged to the true source location and the distances between all robots and the true source location are less than 5 m. Otherwise, it is considered a failure. Additionally, if the search steps reach 30 and the source has not been found, the source search process is terminated as well, indicating an unsuccessful search task. Figure 5 shows the results of single-step search time (SSST) and relative estimation errors (REE) for both algorithms performing a search task using the same simulation parameters as in the example run. It is worth noting that SSST represents the time required to estimate the source position and perform a move at each step, while a REE represents the ratio of the absolute error to actual value.

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Figure 5. PSO-PF vs. PF for multi-robot radiation source search: (A) Relative Estimation Error of PSO-PF; (B) Relative Estimation Error of PF; (C) Single-step search time comparison. EPF: Extended Kalman Particle Filter; PSO-PF: Particle Swarm Optimization-Particle Filter.

Figure 5 presents a clear demonstration of the superior performance of the PSO-PF algorithm in comparison to the EPF algorithm, specifically in terms of the REE and SSST performance indicators. It is important to highlight that the PSO-PF algorithm exhibits a REE below 0.1 for the x-axis coordinate, y-axis coordinate, and activity degree (Ⅰ) of the source term. Conversely, while the EPF algorithm achieves a REE for the source location coordinates below 0.1 at $$ k $$ = 18, the REE for the source activity remains consistently around 0.2. While the PSO-PF algorithm excels in source parameter estimation, it initially incurs additional overhead in terms of single-step search time.

In order to understand the impact of particle numbers on the PSO-PF algorithm, 100 search tasks were also performed using the same parameters as in the example run. Figure 6 shows the search success rate and average search time when varying the number of particles. As the number of particles increases, the search success rate of the free energy strategy increases significantly, climbing from 56% to 98%. At the same time, the average search time increased from 2.75 to 7.04 s. In contrast, the DWA dynamic window obstacle avoidance strategy always maintains a high search success rate as the number of particles increases without significantly affecting the search time. The comparison between the DWA dynamic window obstacle avoidance strategy and the free energy strategy in Figure 6 shows that using the DWA dynamic window obstacle avoidance strategy in the proposed algorithm improves the efficiency and accuracy of multi-robot collaborative search work.

Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

Figure 6. Performance comparison of the PSO-PF algorithm as the number of particles changes (A) The mean search time; (B) The search success rate. PSO-PF: Particle Swarm Optimization-Particle Filter.


In this study, a multi-robot cooperative radioactive source search algorithm based on particle swarm optimization particle filtering is proposed. The algorithm uses the latest observations of each robot to calculate the fitness of the particle swarm, guiding the particle set toward high-likelihood areas. At the same time, the optimized particle set is used in the PF algorithm to achieve accurate source parameter estimation and effectively alleviate the particle degradation problem. In addition, the algorithm combines DWA obstacle avoidance for path planning, allowing the robot to efficiently search for sources in complex mountainous terrains while avoiding obstacles. Experimental results show that the increase in the number of particles increases the search success rate of the algorithm proposed in this study from 90% to 98%, and the average search time increases from 0.97 to 5.38 s. It is worth noting that the radioactive source diffusion and detector model used in this paper is one of the simplified theoretical models. Future research can improve performance by utilizing more realistic models, such as dynamic radiation attenuation coefficients for obstacles, and taking into account factors that affect the intrinsic detection efficiency of the detector model. Lastly, practical experiments using mobile robots equipped with radiation detector modules and incorporating reinforcement learning will be conducted after addressing these issues to bolster system robustness.



The authors would like to thank the editor-in-chief, the associate editor, and the anonymous reviewers for their valuable comments.

Authors' contributions

Made substantial contributions to the research methodology, conceptual, and simulation and wrote and edited the original draft: Luo M, Huo J, Liu M

Performed critical review, comments, and revisions and provided administrative, technical, and material support: Zhou Z, Liu M

Availability of data and materials

Not applicable.

Financial support and sponsorship

This work was supported by the National Natural Science Foundation of China (No. 12205245) and the Natural Science Foundation of Sichuan Province (No. 2023NSFSC1437).

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.


© The Author(s) 2023.


1. IAEA. IAEA releases annual data on illicit trafficking of nuclear and other radioactive material. 2023. Available from: [Last accessed on 18 Dec 2023].

2. Alert in central Mexico after theft of radioactive material. 2023. Available from: [Last accessed on 18 Dec 2023].

3. Howse JW, Ticknor LO, Muske KR. Least squares estimation techniques for position tracking of radioactive sources. Automatica 2001;37:1727-37.

4. Cho HS, Woo TH. Mechanical analysis of flying robot for nuclear safety and security control by radiological monitoring. Ann Nucl Energy 2016;94:138-43.

5. Rao NSV, Shankar M, Chin JC, et al. Identification of low-level point radiation sources using a sensor network. In: 2008 International Conference on Information Processing in Sensor Networks (IPSN 2008); 2008 Apr 22-24; St. Louis, USA. IEEE; 2008. pp. 493-504.

6. Gunatilaka A, Ristic B, Gailis R. On localisation of a radiological point source. In: 2007 Information, Decision and Control; 2007 Feb 12-14; Adelaide, Australia. IEEE; 2007. pp. 236-41.

7. Nemzek RJ, Dreicer JS, Torney DC, Warnock TT. Distributed sensor networks for detection of mobile radioactive sources. IEEE Trans Nucl Sci 2004;51:1693-700.

8. Brewer ET. Autonomous localization of 1/R2 sources using an aerial platform. Virginia Tech; 2009. Available from: [Last accessed on 18 Dec 2023].

9. Morelande MR, Ristic B. Radiological source detection and localisation using Bayesian techniques. IEEE Trans Signal Process 2009;57:4220-31.

10. Jarman KD, Miller EA, Wittman RS, Gesh CJ. Bayesian radiation source localization. Nucl Technol 2011;175:326-34.

11. Huo J, Liu M, Neusypin KA, Liu H, Guo M, Xiao Y. Autonomous search of radioactive sources through mobile robots. Sensors 2020;20:3461.

12. Ling M, Huo J, Moiseev GV, Hu L, Xiao YF. Multi-robot collaborative radioactive source search based on particle fusion and adaptive step size. Ann Nucl Energy 2022;173:109104.

13. Dadgar M, Jafari S, Hamzeh A. A PSO-based multi-robot cooperation method for target searching in unknown environments. Neurocomputing 2016;177:62-74.

14. Couceiro MS, Rocha RP, Ferreira NMF. A novel multi-robot exploration approach based on particle swarm optimization algorithms. In: 2011 IEEE International Symposium on Safety, Security, and Rescue Robotics; 2011 Nov 01-05; Kyoto, Japan. IEEE; 2011. pp. 327-32.

15. Song C, He Y, Ristic B, Lei X. Collaborative infotaxis: searching for a signal-emitting source based on particle filter and Gaussian fitting. Robot Auton Syst 2020;125:103414.

16. Lin HI, Tzeng HJ. Searching a radiological source by a mobile robot. In: 2015 International Conference on Fuzzy Theory and Its Applications (iFUZZY); 2015 Nov 18-20; Yilan, Taiwan. IEEE; 2015. p. 1-5.

17. Gunatilaka A, Ristic B, Gailis R. Radiological source localisation. Defence Science and Technology Organisation Victoria (Australia) Human Protection and Performance Division. 2007. Available from: [Last accessed on 18 Dec 2023].

18. Doucet A, Freitas N, Gordon N. Sequential Monte Carlo methods in practice. vol. 1. In: Information science and statistics. New York: Springer; 2001.

19. Carpenter J, Clifford P, Fearnhead P. Improved particle filter for nonlinear problems. IET Radar Sonar Nav 1999;146:2-7.

20. Kennedy J, Eberhart R. Particle swarm optimization. In: Proceedings of ICNN'95 - International Conference on Neural Networks; 1995 Nov 27 - Dec 01; Perth, Australia. IEEE; 1995. pp. 1942-8.

21. Lu Q, Han QL, Peng D, Choi Y. Decision and event-based fixed-time consensus control for electromagnetic source localization. IEEE Trans Cybern 2022;52:2186-99.

22. Tong G, Fang Z, Xu X. A particle swarm optimized particle filter for nonlinear system state estimation. In: 2006 IEEE International Conference on Evolutionary Computation; 2006 Jul 16-21; Vancouver. IEEE; 2006. pp. 438-42.

23. Zhao J, Li Z. Particle filter based on particle swarm optimization resampling for vision tracking. Expert Syst Appl 2010;37:8910-4.

24. Fox D, Burgard W, Thrun S. The dynamic window approach to collision avoidance. IEEE Robot Autom Mag 1997;4:23-33.

Cite This Article

OAE Style

Luo M, Huo J, Liu M, Zhou Z. Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter. Intell Robot 2023;3(4):685-97.

AMA Style

Luo M, Huo J, Liu M, Zhou Z. Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter. Intelligence & Robotics. 2023; 3(4): 685-97.

Chicago/Turabian Style

Luo, Minghua, Jianwen Huo, Manlu Liu, Zhongbing Zhou. 2023. "Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter" Intelligence & Robotics. 3, no.4: 685-97.

ACS Style

Luo, M.; Huo J.; Liu M.; Zhou Z. Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter. Intell. Robot. 2023, 3, 685-97.

About This Article

Special Issue

This article belongs to the Special Issue Advances in Robotics, AI and Intelligence Control
© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License (, which permits unrestricted use, sharing, adaptation, distribution and reproduction in any medium or format, for any purpose, even commercially, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Data & Comments




Comments must be written in English. Spam, offensive content, impersonation, and private information will not be permitted. If any comment is reported and identified as inappropriate content by OAE staff, the comment will be removed without notice. If you have any queries or need any help, please contact us at

Author's Talk
Download PDF
Cite This Article 6 clicks
Like This Article 7 likes
Share This Article
Scan the QR code for reading!
See Updates
Intelligence & Robotics
ISSN 2770-3541 (Online)
Follow Us


All published articles are preserved here permanently:


All published articles are preserved here permanently: