## Abstract

A swarm robotic system is a system in which multiple robots cooperate to fulfill a macroscopic function. Many swarm robots have been developed for various purposes. This study aims to design swarm robots capable of executing spatially distributed tasks effectively, which can be potentially used for tasks such as search-and-rescue operation and gathering scattered garbage in rooms. We propose a simple decentralized control scheme for swarm robots by extending our previously proposed non-reciprocal-interaction-based model. Each robot has an internal state, called its workload. Each robot first moves randomly to find a task, and when it does, its workload increases, and then it attracts its neighboring robots to ask for their help. We demonstrate, via simulations, that the proposed control scheme enables the robots to effectively execute multiple tasks in parallel under various environments. Fault tolerance of the proposed system is also demonstrated.

## 1 Introduction

Collective behavior emerging from local interaction among individuals is widely found in natural and social systems such as flocking of birds [19, 31] and mammals [5], fish schools [39], bacterial communities [1], and social networks [15]. An interesting aspect of collective behavior is that nontrivial macroscopic functions such as adaptability, scalability, and fault tolerance emerge, although each individual has only trivial functions [29]. While collective behavior has attracted physicists and the mechanism of self-organization has been explored [8, 9, 28], it is also studied from an engineering perspective and many swarm robotic systems have been developed [17, 36].

While previous swarm robots have been developed for many purposes such as aggregation [3, 4, 7, 13, 25, 33, 41], self-assembly [10, 14, 16, 32], collective transport [26], and foraging through task allocation [21, 35, 38, 40], in this study we design a swarm robotic system that can effectively perform spatially distributed tasks. Designing such systems is important because they can potentially be applied to various situations such as search-and-rescue operation, gathering scattered garbage in rooms, planetary surveying, mining, and killing tumor cells in animal bodies [18].

Several previous studies have addressed this issue, some of which assumed that robots can know the locations of tasks before they execute them and scheduled the behaviors of the robots online [2, 6, 11]. However, these studies cannot be applied to cases where robots cannot recognize tasks until they encounter them (for example, in a task of collecting garbage scattered in a room, it is difficult for robots to recognize small garbage before picking it up). Meanwhile, several studies have proposed control schemes that do not require prior information about task locations [3, 12, 20, 30, 33, 34, 41]. However, it is still a challenge to make multiple robots perform tasks fast while minimizing their travel distances under various situations, with reduced computational costs.

Recently, we proposed a non-reciprocal-interaction-based (NRIB) model, and demonstrated, via a simulation, that various patterns emerge through changing parameters [22, 23]. Because of its simplicity, this model has many possible applications in science and engineering, such as understanding the core principle of self-organization [22], elucidating the essential mechanism of the behavior of active matter [37], and designing swarm robotic systems. Thus, the NRIB model is a suitable platform for addressing the above-mentioned issue as well.

This article proposes a decentralized control scheme for an effective execution of spatially distributed tasks by extending the NRIB model. We show, via simulations, that the robots driven by the proposed control scheme effectively execute multiple tasks in parallel in various environments. Moreover, we demonstrate that the proposed control scheme ensures fault tolerance.

The remainder of this article is structured as follows. In Section 2, we briefly review the NRIB model. In Section 3, the problem to be solved in this study is defined. In Section 4, we propose a mathematical model extended from the NRIB model, wherein a decentralized control scheme for the execution of spatially distributed tasks is included. We demonstrate via simulations that the proposed control scheme works well, yet still has a drawback. In Section 5, the control scheme presented in Section 4 is improved to overcome the drawback. The simulation results show that the improved control scheme works well under various conditions. Finally, in Section 6, discussions and conclusions are presented.

## 2 NRIB Model

Let us briefly summarize the NRIB model [22, 23]. Particles, each of which represents a person in a community, exist on a two-dimensional plane, and the position of the ith particle (i = 1, 2, …, N) is denoted by ri. The time evolution of ri is given by
$r˙i=∑j≠ikijRij−1−Rij−2Rˆij,$
(1)
where Rij = rjri, $Rˆ$ij = Rij/|Rij|, and kij denotes a constant that represents to what extent person i prefers person j. The term kij|Rij|−1 in Equation 1 indicates that particle i approaches and repels particle j when kij is positive and negative, respectively. The term −|Rij|−2 in Equation 1 represents the repulsive effect between particles i and j when they are close to each other. When kij = kji, the interaction between the ith and jth particles is described by a potential, and the distance between them tends to converge to kij−1 (if kij > 0). However, because kij is not necessarily equal to kji, that is, the interaction can be non-reciprocal, Equaton 1 is generally a nonequilibrium open system in which both energy and momentum are non-conservative.

The simulation results for several parameter sets of kij can be downloaded from the website link in [22] (http://www.riec.tohoku.ac.jp/%7Etkano/ECAL_Movie1.mp4). It is found that various nontrivial patterns emerge.

## 3 Problem Definition

The problem to be solved in this study is defined in this section. Suppose that N robots are located on a square field of side l. Each robot has an internal state called its workload, and it can move omnidirectionally and detect the relative position and workload of other robots within distance di from itself by using sensors. For simplicity, it is assumed that each robot receives a viscous friction force from the ground. The viscous friction coefficient η is spatially uniform in the field, and is large enough so that the inertia of the robot can be neglected. When two robots collide, they repel each other owing to the exclusive volume effect. The lateral surface of the robot is assumed to be smooth; thus, the friction forces between the robots are neglected. Robots cannot move faster than the maximum velocity vmax. Because this study focuses on extracting the essence of swarm robot control, further details regarding real-world applications (e.g., detailed motor and sensor properties) are not considered.

The tasks to be executed are spatially distributed in the field, and the robots cannot have prior information about their location. When a robot lies in the area where tasks exist, the tasks in the area decrease at a constant rate. At the same time, the workload of the robot increases. All robots have identical properties, and are controlled in a decentralized manner. We will explore a decentralized control scheme for the robots that can execute spatially distributed tasks rapidly with a short travel distance.

## 4 NRIB Swarm Algorithm

Definitions of parameters appearing hereafter are summarized in Tables 1 and 2. The position and workload of robot i (i = 1, 2, …, N) are denoted by ri and wi, respectively. While the space is continuous regarding the robot movement, the amount of tasks in the field is described in a discrete manner. Specifically, the entire field is divided into q × q units, in each of which the amount of tasks z(r) is assumed to be uniform. Each robot can detect the relative positions and workloads of robots that exist in the area within a distance di from itself, denoted by Si. In the present model, di is set to be a constant d0.

Table 1.
Parameter values for conditions 1–3 in the improved NRIB swarm algorithm.
ParameterCondition 1 (nondimensional)Condition 1 (dimensional)
τd 970 970 s
dmin 10.00 10 m
dmax 60 60 m
γ 0.9 0.9
$α˜$′ 0.2 0.2 m · s−1
ParameterCondition 2 (nondimensional)Condition 2 (dimensional)
τd 485 970 s
dmin 5.00 10 m
dmax 30 60 m
γ 0.9 0.9
$α˜$′ 0.2 0.2 m · s−1
ParameterCondition 3 (nondimensional)Condition 3 (dimensional)
τd 323 970 s
dmin 3.33 10 m
dmax 20 60 m
γ 0.9 0.9
$α˜$′ 0.2 0.2 m · s−1
ParameterCondition 1 (nondimensional)Condition 1 (dimensional)
τd 970 970 s
dmin 10.00 10 m
dmax 60 60 m
γ 0.9 0.9
$α˜$′ 0.2 0.2 m · s−1
ParameterCondition 2 (nondimensional)Condition 2 (dimensional)
τd 485 970 s
dmin 5.00 10 m
dmax 30 60 m
γ 0.9 0.9
$α˜$′ 0.2 0.2 m · s−1
ParameterCondition 3 (nondimensional)Condition 3 (dimensional)
τd 323 970 s
dmin 3.33 10 m
dmax 20 60 m
γ 0.9 0.9
$α˜$′ 0.2 0.2 m · s−1
Table 2.
Definition of parameters in the NRIB swarm algorithm and in the improved NRIB swarm algorithm.
VariableDimensionMeaning
N — Number of robots in the field
tmaxstep — Maximum time step
l Length of a side of the square field
q — Number of units for a side of the square field
ri Position of the ith robot
wi — Workload of the ith robot
z(r— Amount of tasks at r
Si — Communication range of the ith robot
di Radius of the communication range (= d0 in Section 4
λ — Coefficient for the magnitude of workload
ϵ s−1 Rate of task execution by a robot
nk — Number of robots that occupy unit k
α kg · m2 · s−2 Coefficient for attractive force
β kg · m3 · s−2 Coefficient for repulsive force
m kg Mass of a robot
η kg · s−1 Viscous friction coefficient between the robot and the ground
f0 kg · m · s−2 Magnitude of the self-driving force of the robot
$fijphys$ kg · m · s−2 Collision force between robots i and j
kphys (depends on μParameter characterizing the collision force
μ — Parameter characterizing the collision force
$α˜$ m2 s−1 = α/η
$β˜$ m3 s−1 = β/η
$f˜$0 m · s−1 = f0/η
$f˜ijphys$ m · s1 = η−1$fijphys$
vmax m · s−1 Maximum velocity of the robot
ni ≡ (cos θi, sin θi)T — Unit vector that points the direction of the self-driving force
td Duration for which the self-driving force is kept constant
tm Maximum value of td
Δt Time step
ET Time until 90% of the tasks are executed
ED Total distance traveled until 90% of the tasks are executed
W — Amount of remaining tasks
D Total distance traveled by all the robots
Parameters for the improved NRIB swarm algorithm
τd Time constant for the change in the communication range
dmin Minimum value of the radius of the communication range
dmax Maximum value of the radius of the communication range
γ — Sensitivity of the communication range to the workload
$α˜$′ m · s−1 Coefficient for attractive force
VariableDimensionMeaning
N — Number of robots in the field
tmaxstep — Maximum time step
l Length of a side of the square field
q — Number of units for a side of the square field
ri Position of the ith robot
wi — Workload of the ith robot
z(r— Amount of tasks at r
Si — Communication range of the ith robot
di Radius of the communication range (= d0 in Section 4
λ — Coefficient for the magnitude of workload
ϵ s−1 Rate of task execution by a robot
nk — Number of robots that occupy unit k
α kg · m2 · s−2 Coefficient for attractive force
β kg · m3 · s−2 Coefficient for repulsive force
m kg Mass of a robot
η kg · s−1 Viscous friction coefficient between the robot and the ground
f0 kg · m · s−2 Magnitude of the self-driving force of the robot
$fijphys$ kg · m · s−2 Collision force between robots i and j
kphys (depends on μParameter characterizing the collision force
μ — Parameter characterizing the collision force
$α˜$ m2 s−1 = α/η
$β˜$ m3 s−1 = β/η
$f˜$0 m · s−1 = f0/η
$f˜ijphys$ m · s1 = η−1$fijphys$
vmax m · s−1 Maximum velocity of the robot
ni ≡ (cos θi, sin θi)T — Unit vector that points the direction of the self-driving force
td Duration for which the self-driving force is kept constant
tm Maximum value of td
Δt Time step
ET Time until 90% of the tasks are executed
ED Total distance traveled until 90% of the tasks are executed
W — Amount of remaining tasks
D Total distance traveled by all the robots
Parameters for the improved NRIB swarm algorithm
τd Time constant for the change in the communication range
dmin Minimum value of the radius of the communication range
dmax Maximum value of the radius of the communication range
γ — Sensitivity of the communication range to the workload
$α˜$′ m · s−1 Coefficient for attractive force

The proposed control scheme is based on “exploration and exploitation” [27]. Specifically, the robots are controlled as follows. Each robot first moves randomly to search for tasks. Once it finds tasks, (i.e., it enters an area with high z(r)) its motion slows down to execute the tasks. Concurrently, its workload increases and its nearby robots detect the increase in the workload and approach it to execute the tasks cooperatively. After the tasks are executed, the workload of the attracted robots decreases. Then, they repel each other and search for other tasks.

Thus, the workload of robot i, wi, evolves according to the following equation:
$τww˙i=λzri−wi,$
(2)
where λ is a positive constant. Equation 2 indicates that wi is given by the first-order delay of λz(ri) with time constant τw. Thus, wi increases when robot i remains in the area with high z(r). Note that an upper limit of wi is introduced; specifically, wi is reset to 1 when it exceeds 1.

It is assumed in the model that the number of tasks decreases at a constant rate in the areas where robots exist. More specifically, z(r) decrements at the rate ϵnk in unit k, where ϵ is a positive constant and nk is the number of robots in unit k. Note that z(r) is reset to zero when it becomes negative. Hence, according to Equation 2, wi decreases as the tasks are executed, that is, z(ri) approaches zero.

The equation of motion for robot i is designed as follows:
$mr¨i+ηr˙i=1−wi∑jϵSigRijRˆij+f0ni+∑jfijphys,$
(3)
where Rij = rjri, $Rˆ$ij = Rij/|Rij|, f0 is a positive constant, and ni is a unit vector. The factor 1 − wi is the mobility of robot i. Specifically, robot i slows down as wi increases, and the driving force completely vanishes when wi = 1. Owing to this factor, robots can continue executing tasks until they are finished. The term f0ni is the self-driving force, which enables random walk of the robot. The unit vector ni is given by ni = (cos θi, sin θi)T. The deflection angle θi takes a uniformly distributed random value in the interval [0, 2π) and is updated after the duration td, where td is a uniformly distributed random number in the interval (0, tm] and is updated when θi is updated. Note that θi is also updated when robot i reaches the edge of the field. When the robot moves out of the field, θi is updated again so that it remains in the field. The force $fijphys$ is the collision force between robots i and j, given by
$fijphys=−kphysmax2r0−Rij0μRˆij,$
(4)
where kphys and μ are positive constants and r0 denotes the radius of the robot.
The function g(|Rij|) determines the interaction between robots. The term g(|Rij|)$Rˆ$ij in Equation 3 denotes attractive and repulsive force when g(|Rij|) is positive and negative, respectively. We designed this function by drawing inspiration from the NRIB model [22, 23]:
$gRij=αwjRij−1−βRij−2,$
(5)
where the first term on the right-hand side indicates that robot i approaches robot j, and its contribution is large when the workload of robot j, wj, is high. Thus, robots tend to aggregate in areas with high z(r). The second term on the right-hand side in Equation 5 means that robot i moves away from robot j when they are close. Note that we designed g(|Rij|) based on the NRIB model, which adopts a power and not an exponential law, because a long-range interaction often helps robots to approach other distant robots.
Equations 3 and 5 are simplified by neglecting the inertia term to become
$r˙i=1−wi∑jϵSiα˜wjRij−1−β˜Rij−2Rˆij+f˜0ni+∑jf˜ijphys,$
(6)
where $α˜$ = α/η, $β˜$ = β/η, $f˜$0 = f0/η, and $f˜ijphys$ = η−1$fijphys$. Note that the absolute value of $r˙$i is reset to vmax when it exceeds vmax.

### 4.1 Simulation Results

Simulations of the proposed model were performed. Areas with high z(r) were randomly distributed under the initial condition. The initial positions of the robots were set to be random. The initial value of the workload wi was set to be 0.2 for all robots.

The simulation program was written in the C++ language. To make the calculation fast without losing accuracy, the time evolution of ri (Equation 6), which often varied rapidly, was solved with the fourth-order Runge-Kutta method, while that of wi (Equation 2), which varied slowly, was solved with the Euler method. Nondimensional parameters were used in the simulation program. The parameters were dimensionalized by multiplying by scaling factors, and three types of scaling factors were examined (Table 3). The parameter values are listed in Table 4. The nondimensional length of a side of the square field was fixed at 40 in the simulation program. This corresponds to 40, 80, and 120 m when a nondimensional length of 1 is scaled to 1, 2, and 3 m, respectively. Meanwhile, the other parameters were determined so that the dimensionalized parameters do not depend on the scaling factors.

Table 3.
Definition of parameters in the BEECLUST algorithm.
VariableDimensionMeaning
p0 — Parameter related to the waiting time
c — Parameter related to the waiting time
ni — Number of robots within the communication range
τbee Duration for the measurement of tr,i, Imin,i, and Imax,i
ncol,i — Number of collisions during the duration τbee
tr,i Average time between subsequent collisions over the duration τbee
$t˜$r,i Average of tr,i over the robots within Si
Imin,i — Minimum value of z(ri) during the duration τbee
Imax,i — Maximum value of z(ri) during the duration τbee
$I˜$min,i — Minimum value of Imin,i for the robots within Si
$I˜$max,i — Maximum value of Imin,i for the robots within Si
Ī — Normalized amount of tasks
δ — Parameter introduced to avoid zero divide (Equation 10
VariableDimensionMeaning
p0 — Parameter related to the waiting time
c — Parameter related to the waiting time
ni — Number of robots within the communication range
τbee Duration for the measurement of tr,i, Imin,i, and Imax,i
ncol,i — Number of collisions during the duration τbee
tr,i Average time between subsequent collisions over the duration τbee
$t˜$r,i Average of tr,i over the robots within Si
Imin,i — Minimum value of z(ri) during the duration τbee
Imax,i — Maximum value of z(ri) during the duration τbee
$I˜$min,i — Minimum value of Imin,i for the robots within Si
$I˜$max,i — Maximum value of Imin,i for the robots within Si
Ī — Normalized amount of tasks
δ — Parameter introduced to avoid zero divide (Equation 10
Table 4.
Scaling factors under conditions 1–3.
ParameterValue
Condition 1Condition 2Condition 3
Nondimensional unit length 1 m 2 m 3 m
Nondimensional unit time 1 s 2 s 3 s
ParameterValue
Condition 1Condition 2Condition 3
Nondimensional unit length 1 m 2 m 3 m
Nondimensional unit time 1 s 2 s 3 s

The simulations were performed 100 times under the same parameters to obtain statistical data. Each trial finished when 90% of the tasks in the field were executed or when the time step reached tmaxstep.

First, simulations were performed with various d0 values under condition 1. Figure 1 shows the time evolution of the remaining tasks W and the total distance D that which all robots traveled, when d0 was 2, 14, and 58 m. The corresponding movies are shown in Movies 1–3 (see the online supplementary material for this article at https://www.mitpressjournals.org/doi/suppl/10.1162/artl_a_00317). When d0 = 2 m, the robots could not approach other robots having high workload, but tended to move independently because of the short communication range. Consequently, the robots had to move randomly for a long time to find tasks; thus, W decreased slowly (Figure 1(a)). When d0 = 58 m, the robots could communicate with all robots in the field. Hence, all robots tended to aggregate into one cluster, so that the tasks were executed cooperatively. However, the robots could not execute multiple tasks in parallel, and hence had to travel long distances; thus, D increased rapidly (Figure 1(c)). In contrast, when d0 = 14 m, the robots formed several clusters to execute multiple tasks in parallel. Once the tasks were executed, the robots were distributed to find other tasks, and formed clusters again when the tasks were found. Thus, W decreased rapidly and D increased slowly, which means that the tasks were executed fast with a short travel distance (Figure 1(b)).

Figure 1.

Time evolution of the remaining tasks W and the total distance for which all robots traveled D when (a) d0 = 2 m, (b) d0 = 14 m, and (c) d0 = 58 m in the NRIB swarm algorithm.

Figure 1.

Time evolution of the remaining tasks W and the total distance for which all robots traveled D when (a) d0 = 2 m, (b) d0 = 14 m, and (c) d0 = 58 m in the NRIB swarm algorithm.

The simulation results were evaluated by the time and the total distance the robots traveled until 90% of the tasks in the field were performed, denoted by ET and ED, respectively. Namely, ET and ED are defined as
$ET=TΔt,$
(7)
$ED=∑k=1T∑i=1NvikΔt,$
(8)
where vi(k) is the velocity vector of robot i at the kth time step and T is the time step when 90% of the tasks in the field are executed, that is, ∫Fz(r)dr decreased by 90% from the initial condition, with ∫F being the spatial integral over the field. Note that T is set to tmaxstep when 90% of the tasks in the field are not executed until the end of the trial.

The result is shown in Figure 2. The median values of ET and ED are the smallest when d0 is 14 and 6 m, respectively. Thus, the robots can execute tasks fast with a short travel distance when d0 is around 10 m.

Figure 2.

Boxplots of ET and ED under condition 1. The results when d0 is varied in the NRIB swarm algorithm and when the improved NRIB swarm algorithm is used are shown. One hundred trials are performed for each d0 value. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars. A statistical test was performed between the NRIB swarm algorithm with d0 =10 m and the improved NRIB swarm algorithm. An asterisk denotes statistically significant difference (p < 0.01). In the statistical analysis, a t-test was adopted when the data is homoscedastic; otherwise a Welch test was adopted.

Figure 2.

Boxplots of ET and ED under condition 1. The results when d0 is varied in the NRIB swarm algorithm and when the improved NRIB swarm algorithm is used are shown. One hundred trials are performed for each d0 value. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars. A statistical test was performed between the NRIB swarm algorithm with d0 =10 m and the improved NRIB swarm algorithm. An asterisk denotes statistically significant difference (p < 0.01). In the statistical analysis, a t-test was adopted when the data is homoscedastic; otherwise a Welch test was adopted.

The result was compared with that of another control scheme proposed previously. Specifically, we compared the result with that of the BEECLUST algorithm [3, 20, 25, 33, 41], a simple decentralized control scheme inspired by the swarming behavior of bees, and evaluated its performance.

The basic concept of the BEECLUST algorithm is as follows:

• 1.

Each robot moves at a constant speed until it collides with an object (another robot, a wall, or an obstacle).

• 2.

When the robot collides with another robot, it stops motion. After the waiting time tbee,i, which is set to become longer as the amount of tasks z(ri) becomes larger, is over, it moves again with a changed direction.

• 3.

When the robot collides with a wall or an obstacle, it immediately changes direction.

The key here is that the waiting time tbee,i is longer as the amount of tasks z(ri) when the collision occurs is larger. Thus, it is expected that the robots stay in the area where the amount of tasks is large, which enables executing tasks effectively.

In the original version of the BEECLUST algorithm [3, 20, 33], robot-to-robot communication is not assumed. That is, no robot knows any information about other robots (such as relative and absolute position). There is no internal memory in the robots. In the extended variant of the BEECLUST algorithm [41], the waiting time tbee,i is adjusted to adapt to the environment by sharing information with other robots. Specifically, information about the robot-to-robot encounter time interval and the maximum and minimum values of z(ri) detected in the recent past is shared with nearby robots. In this study, we performed simulations with the extended variant of the BEECLUST algorithm [41]: The waiting time tbee,i is designed as
$tbee,i=t˜r,ip0I¯i2I¯i2+c,$
(9)
where p0 and c are positive constants, $t˜$r,i denotes the average robot-to-robot encounter time interval, and $I¯$i denotes the normalized amount of tasks. More details about the implementation of the BEECLUST algorithm in this study are provided in the  Appendix.

Figure 3 shows the ET and ED values when p0 is varied. A representative trial when p0 = 5.0 is shown in Movie 4 (see the online supplementary material). For all p0 values, ET and ED are larger than those for the proposed control scheme. According to the statistical analysis, the ET and ED values when d0 = 10 m in the proposed control scheme are significantly smaller than those when p0 = 5 in the BEECLUST algorithm (Welch test, p < 0.01).

Figure 3.

Boxplots of ET and ED when p0 is varied in the BEECLUST algorithm. One hundred trials are performed for each p0 value. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars.

Figure 3.

Boxplots of ET and ED when p0 is varied in the BEECLUST algorithm. One hundred trials are performed for each p0 value. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars.

Next, we performed simulations under conditions 2 and 3, where the spatiotemporal scales are two and three times that under condition 1, respectively (Table 3). The results when d0 is varied are shown in Figure 4. Both ET and ED decrease as d0 increases, and do not vary significantly for large d0. Thus, the optimal value in Condition 1 (d0 ≃ 10 m, Figure 2) is not optimal in conditions 2 and 3. The optimal value of d0 depends on the environment, and the parameters have to be tuned in response to the environment. This is a limitation of the NRIB swarm algorithm.

Figure 4.

Boxplots of ET and ED under (a) condition 1 and (b) condition 2. One hundred trials are performed for each d0 value. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars. A statistical test was performed between the NRIB swarm algorithm with d0 = 10 m and the improved NRIB swarm algorithm. Asterisks denote statistically significant difference (p < 0.01). In the statistical analysis, a t-test was adopted when the data is homoscedastic; otherwise a Welch test was adopted.

Figure 4.

Boxplots of ET and ED under (a) condition 1 and (b) condition 2. One hundred trials are performed for each d0 value. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars. A statistical test was performed between the NRIB swarm algorithm with d0 = 10 m and the improved NRIB swarm algorithm. Asterisks denote statistically significant difference (p < 0.01). In the statistical analysis, a t-test was adopted when the data is homoscedastic; otherwise a Welch test was adopted.

## 5 Improved NRIB Swarm Algorithm

To solve the problem of the NRIB swarm algorithm, an improved NRIB swarm algorithm is proposed, where the radius of the communication range di changes dynamically. The basic idea is as follows. When there does not exist a robot having high workload within the communication range, the robot needs to enlarge the communication range to find robots that need help. In contrast, when there exist many robots having high workload within the communication range, the robot will become confused about which direction to move in. In such a case, it is desirable that the robot approach the nearest robot with high workload so as not to move a long distance; hence, the communication range needs to be reduced.

Based on the above idea, the time evolution of di is designed as
$τrd˙i=dmax−dmax−dmintanhγ∑jϵSiwj−di,$
(10)
where γ and τr are positive constants. The radius of the communication range di approaches dmax and dmin as the summation of the workload of the robots within the communication range increases and decreases, respectively. Thus, it is expected that the communication range will be auto-tuned in response to the situation.
Moreover, the algorithm for the movement of the robots is slightly modified from Equation 6. When the communication range di is large, robot i is required to approach distant robots with high workload to execute tasks cooperatively. Thus, the attraction term ($α˜$wi|Rij|−1 in Equation 6) should increase as di increases. Hence, Equation 6 is modified as
$r˙i=1−wi∑jϵSiα˜′diwjRij−1−β˜Rij−2Rˆij+f˜0ni+∑jf˜ijphys,$
(11)
where $α˜$′ is a positive constant.

### 5.1 Simulation Results

Simulations were performed using the improved NRIB swarm algorithm. The time evolution of di (Equation 10) was solved using the Euler method, because it varied slowly. The values of parameters specific to the improved NRIB swarm algorithm are listed in Table 5, and the other parameter values are the same as those of the NRIB swarm algorithm (Table 4). The results for conditions 1–3 are shown in Figures 2, 4(a), and 4(b) and Movies 5–7 (see the online supplementary material). It is found that both ET and ED for the improved NRIB swarm algorithm are smaller than or almost equal to the optimal values for the NRIB swarm algorithm. According to statistical analysis, ET and ED for the improved NRIB swarm algorithm are significantly smaller than those for the NRIB swarm algorithm with d0 = 10 m (Welch test or t-test, p < 0.01), except for ED for condition 1, wherein no significant difference was found. This was achieved by auto-tuning of the communication range.

Table 5.
Parameter values under conditions 1–3 in the NRIB swarm algorithm.
ParameterValue
Condition 1 (non-dimensional)Condition 1 (dimensional)Condition 2 (non-dimensional)Condition 2 (dimensional)Condition 3 (non-dimensional)Condition 3 (dimensional)
N 30 30 30 30 30 30
tmaxstep 2000000 2000000 2000000 2000000 2000000 2000000
r0 0.500 0.5 m 0.250 0.5 m 0.167 0.5 m
l 40 40 m 40 80 m 40 120 m
q 200 200 200 200 200 200
τw 1.000 1 s 0.500 1 s 0.333 1 s
λ
ϵ 0.1 0.1 s−1 0.2 0.1 s−1 0.3 0.1 s−1
$α˜$ 2.000 2 m2 s−1 1.000 2 m2 s−1 0.667 2 m2 s−1
$β˜$ 2.000 2 m3 s−1 0.500 2 m3 s−1 0.222 2 m3 s−1
$f˜$0 0.5 0.5 m · s−1 0.5 0.5 m · s−1 0.5 0.5 m · s−1
$k˜$phys(≡ kphys/η1000 1000 m−4 s−1 32000 1000 m−4 s−1 243000 1000 m−4 s−1
μ
vmax 1 m · s−1 1 m · s−1 1 m · s−1
tm 300 300 s 150 300 s 100 300 s
Δt 0.001000 0.001 s 0.000500 0.001 s 0.000333 0.001 s
ParameterValue
Condition 1 (non-dimensional)Condition 1 (dimensional)Condition 2 (non-dimensional)Condition 2 (dimensional)Condition 3 (non-dimensional)Condition 3 (dimensional)
N 30 30 30 30 30 30
tmaxstep 2000000 2000000 2000000 2000000 2000000 2000000
r0 0.500 0.5 m 0.250 0.5 m 0.167 0.5 m
l 40 40 m 40 80 m 40 120 m
q 200 200 200 200 200 200
τw 1.000 1 s 0.500 1 s 0.333 1 s
λ
ϵ 0.1 0.1 s−1 0.2 0.1 s−1 0.3 0.1 s−1
$α˜$ 2.000 2 m2 s−1 1.000 2 m2 s−1 0.667 2 m2 s−1
$β˜$ 2.000 2 m3 s−1 0.500 2 m3 s−1 0.222 2 m3 s−1
$f˜$0 0.5 0.5 m · s−1 0.5 0.5 m · s−1 0.5 0.5 m · s−1
$k˜$phys(≡ kphys/η1000 1000 m−4 s−1 32000 1000 m−4 s−1 243000 1000 m−4 s−1
μ
vmax 1 m · s−1 1 m · s−1 1 m · s−1
tm 300 300 s 150 300 s 100 300 s
Δt 0.001000 0.001 s 0.000500 0.001 s 0.000333 0.001 s

We also performed simulations under condition 1 with the following cases: (i) tasks appeared randomly in the field during the simulation, (ii) some robots suddenly stopped moving during the simulation, and (iii) obstacles existed in the field. Note that maximal tmaxstep was set to be 100000 in case (ii), while it was 200000 in the other cases. Representative trials for cases (i)–(iii) are shown in Movies 8–10 in the supplementary material. We found that the robots could execute the tasks for all cases.

The results for cases (i)–(iii) with the improved NRIB swarm algorithm were compared with those for the NRIB swarm algorithm and the BEECLUST algorithm. The radius of the communication range d0 was set to 10 m in the NRIB swarm algorithm, and p0 was set to 5.0 in the BEECLUST algorithm. The performance was evaluated by ET and ED in cases (ii) and (iii); however, ET could not be used in case (i) because of the appearance of tasks. Hence, the performance in case (i) was evaluated by ED and ER, the latter of which is defined as the ratio of ∫Fz(r)dr at t = tmaxstep to that at t = 0.

The results are shown in Figure 5. In case (i), ER and ED for the BEECLUST algorithm are significantly larger than for the other algorithms (Welch test, p < 0.01). ER for the improved NRIB swarm algorithm is significantly smaller than that for the NRIB swarm algorithm (Welch test, p < 0.01). However, ED for the improved NRIB swarm algorithm is significantly larger than that for the NRIB swarm algorithm (Welch test, p < 0.01). This result was obtained because robots moved toward distant robots to execute tasks fast and cooperatively (Movie 8). In cases (ii) and (iii), ET and ED for the improved NRIB swarm algorithm are significantly smaller than those for the BEECLUST algorithm and the NRIB swarm algorithm (Welch test or t-test, p < 0.01 or p < 0.05).

Figure 5.

Boxplots of ER, ET, and ED when (i) tasks appear randomly, (ii) some robots stop moving, and (iii) obstacles exist. Comparisons among the BEECLUST algorithm (p0 = 5), the NRIB swarm algorithm (d0 = 10 m), and the improved NRIB swarm algorithm are shown. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars. The ET value for the BEECLUST algorithm in (ii) is not shown, because 90% of the tasks were not executed until the maximum time step tmaxstep. Asterisks denote statistically significant difference. In the statistical analysis, the t-test was adopted when the data is homoscedastic; otherwise the Welch test was adopted.

Figure 5.

Boxplots of ER, ET, and ED when (i) tasks appear randomly, (ii) some robots stop moving, and (iii) obstacles exist. Comparisons among the BEECLUST algorithm (p0 = 5), the NRIB swarm algorithm (d0 = 10 m), and the improved NRIB swarm algorithm are shown. The median and quartile are shown by boxes, while the maximum and minimum values are shown by bars. The ET value for the BEECLUST algorithm in (ii) is not shown, because 90% of the tasks were not executed until the maximum time step tmaxstep. Asterisks denote statistically significant difference. In the statistical analysis, the t-test was adopted when the data is homoscedastic; otherwise the Welch test was adopted.

## 6 Discussion and Conclusion

We have proposed, by drawing inspiration from the NRIB model proposed previously [22, 23], simple decentralized control schemes for swarm robots that can execute spatially distributed tasks in parallel. We proposed two control schemes (the NRIB swarm algorithm and the improved NRIB swarm algorithm). In the NRIB swarm algorithm (Section 4), each robot moves randomly to find tasks; once it finds the tasks, it attracts robots within its communication range to execute them cooperatively. We performed simulations under various conditions and found that the communication range d0 is an important parameter. For small d0, the robots move almost independently and need to travel for a long time to find tasks, while for large d0 they move long distances to aggregate (Figure 1). Thus, d0 has an optimal value; however, it depends on the length scale of the field (conditions 1–3, Figures 2 and 4). To solve this problem, we proposed an improved NRIB swarm algorithm, wherein the communication range is auto-tuned on the basis of the workload of robots within the communication range (Section 5). This auto-tuning mechanism is reasonable because a robot searches a large area when robots that need help do not exist in the vicinity. We demonstrated, via simulations, that the improved NRIB swarm algorithm is applicable to various situations without changing any parameter (Figures 2, 4, and 5).

We compared the proposed control schemes with the extended variant of the BEECLUST algorithm [41]. The BEECLUST algorithm is simpler than the proposed control schemes and does not need to detect relative positions of other robots: Robots only need to share information about the robot-to-robot encounter time interval and the maximum and minimum values of z(ri) detected in the recent past. Whereas the BEECLUST algorithm has this advantage, the proposed control scheme enables the robots to execute tasks fast with a short travel distance compared with the BEECLUST algorithm (Figures 3 and 5). This result was obtained because robots performed tasks cooperatively only when they collided with each other by coincidence in the BEECLUST algorithm, whereas the proposed control schemes enabled cooperative task execution through interaction between robots.

Several algorithms besides the BEECLUST algorithm have been proposed thus far [12, 30, 34]. Although we have not quantitatively compared the proposed control schemes with these studies, the idea of the variable communication range proposed in the improved NRIB swarm algorithm (Equation 10) is novel, and is expected to be advantageous. In fact, we showed that the performance improved by auto-tuning of the communication range (Figures 2 and 4).

We believe that the proposed control scheme can be used for various practical applications. Hardware realization is a first step toward this. We believe that it is possible in the near future, because we have already succeeded in developing hardware for the NRIB model wherein sensors that can detect the relative position of nearby robots and actuators that enable omnidirectional movement are incorporated [24].

The proposed control schemes still have limitations. First, they may not be applicable to the cases where the length scale of tasks is comparable with or smaller than the robot diameter. This is because some robots with high workload occupy the area where a task exists, and thus, robots attracted by them cannot enter there. Second, the performance of the proposed control scheme is not guaranteed when the ground friction is inhomogeneous. For example, when there exists an area with extremely high friction, robots need to avoid entering it; otherwise, they can get stuck; however, the present control schemes cannot cope with such a situation. Further studies are needed to solve these problems.

## Acknowledgment

The authors thank Professor Ken Sugawara of Tohoku Gakuin University and Naoki Matsui of Tohoku University for their helpful suggestions.

## References

1
Asally
,
M.
,
Kittisopikul
,
M.
,
Rué
,
P.
,
Du
,
Y.
,
Hu
,
Z.
,
Cağatay
,
T.
,
Robinson
,
A. B.
,
Lu
,
H.
,
Garcia-Ojalvo
,
J.
, &
Süel
,
G. M.
(
2012
).
Localized cell death focuses mechanical forces during 3D patterning in a biofilm
.
Proceedings of the National Academy of Sciences of the USA
,
109
(
46
),
18891
18896
. https://doi.org/10.1073/pnas.1212429109.
2
Aşık
,
O.
, &
Akın
,
H. L.
(
2017
).
Effective multi-robot spatial task allocation using model approximations
. In
S.
Behnke
,
R.
Sheh
,
S.
Sarıel
, &
D.
Lee
(Eds.),
RoboCup 2016: Robot World Cup XX
(pp.
243
255
).
Cham, Switzerland
:
Springer
. https://doi.org/10.1007/978-3-319-68792-6_20.
3
Bodi
,
M.
,
Thenius
,
R.
,
Szopek
,
M.
,
Schmickl
,
T.
, &
Crailsheim
,
K.
(
2012
).
Interaction of robot swarms using the honeybee-inspired control algorithm BEECLUST
.
Mathematical and Computer Modelling of Dynamical Systems
,
18
(
1
),
87
100
. https://doi.org/10.1080/13873954.2011.601420.
4
Cambier
,
N.
,
Fŕemont
,
V.
, &
Ferrante
,
E.
(
2017
).
Group-size regulation in self-organised aggregation through the naming game
. In
Proceedings of International Symposium on Swarm Behavior and Bio-Inspired Robotics (SWARM 2017)
(pp.
365
372
).
5
Chen
,
Y.
, &
Kolokolnikov
,
T.
(
2014
).
A minimal model of predator-swarm interactions
,
Journal of the Royal Society Interface
,
11
,
20131208
. https://doi.org/10.1098/rsif.2013.1208.
6
Claes
,
D.
,
Robbel
,
P.
,
Oliehoek
,
F. A.
,
Tuyls
,
K.
,
Hennes
,
D.
, &
van der Hoek
,
W.
(
2015
).
Effective approximations for multi-robot coordination in spatially distributed tasks
. In
B.
Bordini
,
M. S. V.
Elkind
,
R. G.
Weiss
, &
G. T.
Yolum
(Eds.),
Proceedings of the 2015 International Conference on Autonomous Agents and Multiagent Systems
(pp.
881
890
).
Richland, SC
:
International Foundation for Autonomous Agents and Multiagent Systems (IFAAMAS)
.
7
Correll
,
N.
, &
Martinoli
,
A.
(
2011
).
Modeling and designing self-organized aggregation in a swarm of miniature robots
.
International Journal of Robotic Research
,
30
(
5
),
615
626
. https://doi.org/10.1177/0278364911403017.
8
Czirók
,
A.
,
Stanley
,
H. E.
, &
Vicsek
,
T.
(
1997
).
Spontaneously ordered motion of self-propelled particles
.
Journal of Physics A: Mathematical and General
,
30
(
5
),
1375
. https://doi.org/10.1088/0305-4470/30/5/009.
9
Czirók
,
A.
, &
Vicsek
,
T.
(
2000
).
Collective behavior of interacting self-propelled particles
.
Physica A: Statistical Mechanics and its Applications
,
281
(
1–4
),
17
29
. https://doi.org/10.1016/S0378-4371(00)00013-3.
10
Dorigo
,
M.
,
Tuci
,
E.
,
Trianni
,
V.
Groß
,
R.
,
Nouyan
,
R.
,
Ampatzis
,
C.
,
Labella
,
T. H.
,
,
R.
,
Bonani
,
M.
, &
F.
(
2006
).
SWARMBOT: Design and implementation of colonies of self-assembling robots
. In
G. Y.
Yen
&
D. B.
Fogel
(Eds.),
Computational intelligence: Principles and practice
(pp.
103
135
).
Los Alamitos, CA
:
IEEE Press
.
11
Ducatelle
,
F.
,
Förster
,
A.
,
Di Caro
,
G. A.
, &
Gambardella
,
L. M.
(
2009
).
New task allocation methods for robotic swarms
. In
Proceedings of 9th IEEE/RAS Conference on Autonomous Robot Systems and Competitions
.
12
Dasgupta
,
P.
(
2011
).
Multi-robot task allocation for performing cooperative foraging tasks in an initially unknown environment
. In
L. C.
Jain
et al
(Eds.),
Innovations in defence support systems—2, SCI 338
(pp.
5
20
).
Berlin, Heidelberg
:
Springer Verlag
. https://doi.org/10.1007/978-3-642-17764-4.
13
Firat
,
Z.
,
Ferrante
,
E.
,
Gilet
,
Y.
, &
Tuci
,
E.
(
2019
).
On self-organised aggregation dynamics in swarms of robots with informed robots
.
ArXiv:1903.03841
.
14
Garnier
,
S.
,
Murphy
,
T.
,
Lutz
,
M.
,
Hurme
,
E.
,
Leblanc
,
S.
, &
Couzin
,
I. D.
(
2013
).
Stability and responsiveness in a self-organized living architecture
.
PLOS Computational Biology
,
9
(
3
),
e1002984
. https://doi.org/10.1371/journal.pcbi.1002984.
15
González
,
M. C.
,
Lind
,
P. G.
, &
Herrman
,
H. J.
(
2006
).
System of mobile agents to model social networks
.
Physical Review Letters
,
96
,
088702
. https://doi.org/10.1103/physrevlett.96.088702.
16
Groß
,
R.
,
Bonani
,
M.
,
,
F.
, &
Dorigo
,
M.
(
2006
).
Autonomous self-assembly in swarm-bots
.
IEEE Transactions on Robotics
,
22
(
6
),
1115
1130
. https://doi.org/10.1109/TRO.2006.882919.
17
Hamann
,
H.
(
2018
).
Swarm robotics: A formal approach
(1st ed.).
New York
:
Springer
.
18
Hauert
,
S.
, &
Bhatia
,
S. N.
(
2014
).
Mechanisms of cooperation in cancer nanomedicine: Towards systems nanotechnology
.
Trends in Biotechnology
,
32
,
448
455
. https://doi.org/10.1016/j.tibtech.2014.06.010.
19
Hayakawa
,
Y.
(
2010
).
Spatiotemporal dynamics of skeins of wild geese
.
Europhysics Letters
,
89
,
48004
. https://doi.org/10.1209/0295-5075/89/48004.
20
Hereford
,
J.
(
2011
).
Analysis of BEECLUST swarm algorithm
. In
Proceedings of IEEE Symposium on Swarm Intelligence (SIS)
(pp.
1
7
).
New York
:
IEEE
. https://doi.org/10.1109/SIS.2011.5952587.
21
Jamshidpey
,
A.
, &
Afsharchi
,
M.
(
2015
).
Task allocation in robotic swarms: Explicit communication based approaches
. In
D.
Barbosa
&
E.
Milios
(Eds.),
Proceedings of Canadian Conference on Artificial Intelligence
(pp.
59
67
).
New York
:
Springer
. https://doi.org/10.1007/978-3-319-18356-5_6.
22
Kano
,
T.
,
Osuka
,
K.
,
Kawakatsu
,
T.
,
Matsui
,
N.
, &
Ishiguro
,
A.
(
2017
).
A minimal model of collective behaviour based on non-reciprocal Interactions
. In
C.
Knibbe
,
G.
Beslon
,
D. P.
Parsons
,
D.
Misevic
,
J.
Rouzaud-Cornabas
,
N.
Bredèche
,
S.
Hassas
,
O.
Simonin
, &
H.
Soula
(Eds.),
Proceedings of ECAL 2017 the 14th European Conference on Artificial Life
(pp.
237
244
).
Cambridge, MA
:
MIT Press
. https://doi.org/10.7551/ecal_a_041.
23
Kano
,
T.
,
Osuka
,
K.
,
Kawakatsu
,
T.
, &
Ishiguro
,
A.
(
2017
).
Mathematical analysis for non-reciprocal-interaction-based model of collective behavior
.
Journal of the Physical Society of Japan
,
86
,
124004
. https://doi.org/10.7566/JPSJ.86.124004.
24
Kano
,
T.
,
Matsui
,
N.
,
Naito
,
E.
,
Aoshima
,
T.
, &
Ishiguro
,
A.
(
2018
).
Swarm robots inspired by friendship formation process
.
arXiv:1808.03812
.
25
Kernbach
,
S.
,
Häbe
,
D.
,
Kernbach
,
O.
,
Thenius
,
R.
,
,
G.
,
Kimura
,
T.
, &
Schmickl
,
T.
(
2012
).
Adaptive collective decision-making in limited robot swarms without communication
.
The International Journal of Robotics Research
,
32
(
1
),
35
55
. https://doi.org/10.1177/0278364912468636.
26
Kube
C. R.
, &
Bonabeau
,
E.
(
2000
).
Cooperative transport by ants and robots
.
Robotics and Autonomous Systems
,
30
,
85
101
. https://doi.org/10.1016/S0921-8890(99)00066-4.
27
March
,
J. G.
(
1991
).
Exploration and exploitation in organizational learning
.
Organization Science
,
2
,
71
87
. https://doi.org/10.1287/orsc.2.1.71.
28
Mishra
,
S.
,
,
A.
, &
Marchetti
,
M. C.
(
2010
).
Fluctuations and pattern formation in self-propelled particles
.
Physical Review E
,
81
(
6
),
061916
. https://doi.org/10.1103/PhysRevE.81.061916.
29
Petersen
,
K. H.
,
Napp
,
N.
,
Stuart-Smith
,
R.
,
Rus
,
D.
, &
Kovac
,
M.
(
2019
).
A review of collective robotic construction
.
Science Robotics
,
4
,
eaau8479
. https://doi.org/10.1126/scirobotics.aau8479.
30
Reif
,
J. H.
, &
Wang
,
H.
(
1999
).
Social potential fields: A distributed behavioral control for autonomous robots
.
Robotics and Autonomous Systems
,
27
,
171
194
. https://doi.org/10.1016/S0921-8890(99)00004-4.
31
Reynolds
,
C. W.
(
1987
).
Flocks, herds, and schools: A distributed behavioral model
.
Computer Graphics
,
21
,
25
34
. https://doi.org/10.1145/37402.37406.
32
Rubenstein
,
M.
,
Cornejo
,
A.
, &
Napgal
,
R.
(
2014
).
Programmable self-assembly in a thousand-robot swarm
.
Science
,
345
(
6198
),
795
799
. https://doi.org/10.1126/science.1254295.
33
Schmickl
,
T.
, &
Hamann
,
H.
(
2011
).
BEECLUST: A swarm algorithm derived from honeybees
. In
X. S.
Yang
et al
(Eds.),
Bio-inspired computing and communication networks
.
Boca Raton, FL
:
CRC Press
.
34
Simonin
,
O.
, &
Ferber
,
J.
(
2000
).
Modeling self satisfaction and altruism to handle action selection and reactive cooperation
. In
J.
Meyer
,
A.
Berthoz
,
D.
Floreano
,
H.
Roitblat
, &
S. W.
Wilson
(Eds.),
Proceedings of 6th International Conference on the Simulation of Adaptive Behavior
(pp.
314
323
).
Cambridge, MA
:
MIT Press
.
35
Sugawara
,
K.
, &
Sano
,
M.
(
1997
).
Cooperative acceleration of task performance: Foraging behavior of interacting multi-robots system
.
Physica D
,
100
,
343
354
. https://doi.org/10.1016/S0167-2789(96)00195-9.
36
Tan
,
Y.
, &
Zheng
,
Z.
(
2013
).
.
Defence Technology
,
9
,
18
39
. https://doi.org/10.1016/j.dt.2013.03.001.
37
Tanaka
,
S.
,
Nakata
,
S.
, &
Kano
,
T.
(
2017
).
Dynamic ordering in a swarm of floating droplets driven by solutal Marangoni effect
.
Journal of the Physical Society of Japan
,
86
,
101004
. https://doi.org/10.7566/JPSJ.86.101004.
38
Theraulaz
,
G.
,
Bonabeau
,
E.
, &
Deneubourg
,
J.-L.
(
1998
).
Response threshold reinforcements and division of labour in insect societies
.
Proceedings of the Royal Society London: Series B
,
265
,
327
332
. https://doi.org/10.1098/rspb.1998.0299.
39
Vabø
,
R.
, &
,
L.
(
1997
).
An individual based model of fish school reactions: Predicting antipredator behaviour as observed in nature
.
Fisheries Oceanography
,
6
(
3
),
155
171
. https://doi.org/10.1046/j.1365-2419.1997.00037.x.
40
Van Essche
,
S.
,
Ferrante
E.
,
Turgut
,
A. E.
,
Van Lon
,
R.
,
Holvoet
,
T.
, &
Wenseleers
,
T.
(
2015
).
Environmental factors promoting the evolution of recruitment strategies in swarms of foraging robots
. In
Proceedings of International Symposium on Swarm Behavior and Bio-Inspired Robotics (SWARM 2015)
(pp.
389
396
).
41
Wahby
,
M.
,
Petzold
,
J.
,
Eschke
,
C.
,
Schmickl
,
T.
, &
Hamann
,
H.
(
2019
).
Collective change detection: Adaptivity to dynamic swarm densities and light conditions in robot swarms
. In
H.
Fellermann
,
J.
Bacardit
,
Á. G.
Moreno
, &
R.
Füchslin
(Eds.),
Proceedings on Artificial Life Conference
(pp.
642
649
).
Cambridge, MA
:
MIT Press
. https://doi.org/10.1162/isal_a_00233.

### Appendix

Here we explain how the extended variant of the BEECLUST algorithm [41] was implemented in our simulator. A flowchart of the algorithm is shown in Figure 6.

Figure 6.

Flowchart of the the extended variant of the BEECLUST algorithm [41] implemented in our simulator.

Figure 6.

Flowchart of the the extended variant of the BEECLUST algorithm [41] implemented in our simulator.

Basically, each robot moves with the velocity vector vmaxni, where ni = (cos θi, sin θi)T. The deflection angle θi is given by a random number in the interval [0, 2π) and remains unchanged until the robot collides with an object (another robot, a wall, or an obstacle).

Each robot always measures tr,i, Imin,i, and Imax,i defined below and shares them with other robots within its communication range.
• tr,i: Average time between subsequent collisions over the duration τbee in the recent past. It is defined by tr,i = τbee/ncol,i, where ncol,i denotes the number of collisions during the duration τbee.

• Imin,i: Minimum value of z(ri) during the duration τbee in the recent past.

• Imax,i: Maximum value of z(ri) during the duration τbee in the recent past.

By sharing this information, each robot obtains $t˜$r,i, $I˜$min,i, and $I˜$max,i given by
$t˜r,i=1ni∑jϵSitr,j,I˜min,i=minjϵSiImin,j,I˜max,i=maxjϵSiImax,jI˜min,i+δ,$
(12)
where ni denotes the number of robots within the communication range, and δ is a positive constant that was introduced to avoid zero divide in Equation 13 described below. Then, the normalized amount of tasks $I˜$i is calculated as
$I¯i=zri−I˜min,iI˜max,i−I˜min,i$
(13)
Using $t˜$r,i and $I¯$i, the waiting time tbee,i is given by Equation 9 in the main text.

When the robot collides with another robot, it stops motion for the duration tbee and then moves with the velocity vector vmaxni again, updating the θi value randomly. When the robot collides with a wall or an obstacle, it immediately changes the direction of motion by updating the θi value randomly.