## 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

*i*th particle (

*i*= 1, 2, …,

*N*) is denoted by

**r**

_{i}. The time evolution of

**r**

_{i}is given by

**R**

_{ij}=

**r**

_{j}−

**r**

_{i}, $R\u02c6$

_{ij}=

**R**

_{ij}/|

**R**

_{ij}|, and

*k*

_{ij}denotes a constant that represents to what extent person

*i*prefers person

*j*. The term

*k*

_{ij}|

**R**

*ij*|

^{−1}in Equation 1 indicates that particle

*i*approaches and repels particle

*j*when

*k*

_{ij}is positive and negative, respectively. The term −|

**R**

_{ij}|

^{−2}in Equation 1 represents the repulsive effect between particles

*i*and

*j*when they are close to each other. When

*k*

_{ij}=

*k*

_{ji}, the interaction between the

*i*th and

*j*th particles is described by a potential, and the distance between them tends to converge to

*k*

_{ij}

^{−1}(if

*k*

_{ij}> 0). However, because

*k*

_{ij}is not necessarily equal to

*k*

_{ji}, 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 *k*_{ij} 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 *d*_{i} 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 *v*_{max}. 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 **r**_{i} and *w*_{i}, 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 *d*_{i} from itself, denoted by *S*_{i}. In the present model, *d*_{i} is set to be a constant *d*_{0}.

Parameter . | Condition 1 (nondimensional) . | Condition 1 (dimensional) . |
---|---|---|

τ_{d} | 970 | 970 s |

d_{min} | 10.00 | 10 m |

d_{max} | 60 | 60 m |

γ | 0.9 | 0.9 |

$\alpha \u02dc$′ | 0.2 | 0.2 m · s^{−1} |

Parameter . | Condition 2 (nondimensional) . | Condition 2 (dimensional) . |

τ_{d} | 485 | 970 s |

d_{min} | 5.00 | 10 m |

d_{max} | 30 | 60 m |

γ | 0.9 | 0.9 |

$\alpha \u02dc$′ | 0.2 | 0.2 m · s^{−1} |

Parameter . | Condition 3 (nondimensional) . | Condition 3 (dimensional) . |

τ_{d} | 323 | 970 s |

d_{min} | 3.33 | 10 m |

d_{max} | 20 | 60 m |

γ | 0.9 | 0.9 |

$\alpha \u02dc$′ | 0.2 | 0.2 m · s^{−1} |

Parameter . | Condition 1 (nondimensional) . | Condition 1 (dimensional) . |
---|---|---|

τ_{d} | 970 | 970 s |

d_{min} | 10.00 | 10 m |

d_{max} | 60 | 60 m |

γ | 0.9 | 0.9 |

$\alpha \u02dc$′ | 0.2 | 0.2 m · s^{−1} |

Parameter . | Condition 2 (nondimensional) . | Condition 2 (dimensional) . |

τ_{d} | 485 | 970 s |

d_{min} | 5.00 | 10 m |

d_{max} | 30 | 60 m |

γ | 0.9 | 0.9 |

$\alpha \u02dc$′ | 0.2 | 0.2 m · s^{−1} |

Parameter . | Condition 3 (nondimensional) . | Condition 3 (dimensional) . |

τ_{d} | 323 | 970 s |

d_{min} | 3.33 | 10 m |

d_{max} | 20 | 60 m |

γ | 0.9 | 0.9 |

$\alpha \u02dc$′ | 0.2 | 0.2 m · s^{−1} |

Variable . | Dimension . | Meaning . |
---|---|---|

N | — | Number of robots in the field |

t_{maxstep} | — | Maximum time step |

r_{0} | m | Radius of the robot |

l | m | Length of a side of the square field |

q | — | Number of units for a side of the square field |

r_{i} | m | Position of the ith robot |

w_{i} | — | Workload of the ith robot |

z(r) | — | Amount of tasks at r |

S_{i} | — | Communication range of the ith robot |

d_{i} | m | Radius of the communication range (= d_{0} in Section 4) |

τ_{w} | s | Time constant for workload |

λ | — | Coefficient for the magnitude of workload |

ϵ | s^{−1} | Rate of task execution by a robot |

n_{k} | — | Number of robots that occupy unit k |

α | kg · m^{2} · s^{−2} | Coefficient for attractive force |

β | kg · m^{3} · s^{−2} | Coefficient for repulsive force |

m | kg | Mass of a robot |

η | kg · s^{−1} | Viscous friction coefficient between the robot and the ground |

f_{0} | 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 |

k^{phys} | (depends on μ) | Parameter characterizing the collision force |

μ | — | Parameter characterizing the collision force |

$\alpha \u02dc$ | m^{2} s^{−1} | = α/η |

$\beta \u02dc$ | m^{3} s^{−1} | = β/η |

$f\u02dc$_{0} | m · s^{−1} | = f_{0}/η |

$f\u02dcijphys$ | m · s^{−1} | = η^{−1}$fijphys$ |

v_{max} | m · s^{−1} | Maximum velocity of the robot |

n_{i} ≡ (cos θ_{i}, sin θ_{i})^{T} | — | Unit vector that points the direction of the self-driving force |

t_{d} | s | Duration for which the self-driving force is kept constant |

t_{m} | s | Maximum value of t_{d} |

Δt | s | Time step |

E_{T} | s | Time until 90% of the tasks are executed |

E_{D} | m | Total distance traveled until 90% of the tasks are executed |

W | — | Amount of remaining tasks |

D | m | Total distance traveled by all the robots |

Parameters for the improved NRIB swarm algorithm | ||

τ_{d} | s | Time constant for the change in the communication range |

d_{min} | m | Minimum value of the radius of the communication range |

d_{max} | m | Maximum value of the radius of the communication range |

γ | — | Sensitivity of the communication range to the workload |

$\alpha \u02dc$′ | m · s^{−1} | Coefficient for attractive force |

Variable . | Dimension . | Meaning . |
---|---|---|

N | — | Number of robots in the field |

t_{maxstep} | — | Maximum time step |

r_{0} | m | Radius of the robot |

l | m | Length of a side of the square field |

q | — | Number of units for a side of the square field |

r_{i} | m | Position of the ith robot |

w_{i} | — | Workload of the ith robot |

z(r) | — | Amount of tasks at r |

S_{i} | — | Communication range of the ith robot |

d_{i} | m | Radius of the communication range (= d_{0} in Section 4) |

τ_{w} | s | Time constant for workload |

λ | — | Coefficient for the magnitude of workload |

ϵ | s^{−1} | Rate of task execution by a robot |

n_{k} | — | Number of robots that occupy unit k |

α | kg · m^{2} · s^{−2} | Coefficient for attractive force |

β | kg · m^{3} · s^{−2} | Coefficient for repulsive force |

m | kg | Mass of a robot |

η | kg · s^{−1} | Viscous friction coefficient between the robot and the ground |

f_{0} | 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 |

k^{phys} | (depends on μ) | Parameter characterizing the collision force |

μ | — | Parameter characterizing the collision force |

$\alpha \u02dc$ | m^{2} s^{−1} | = α/η |

$\beta \u02dc$ | m^{3} s^{−1} | = β/η |

$f\u02dc$_{0} | m · s^{−1} | = f_{0}/η |

$f\u02dcijphys$ | m · s^{−1} | = η^{−1}$fijphys$ |

v_{max} | m · s^{−1} | Maximum velocity of the robot |

n_{i} ≡ (cos θ_{i}, sin θ_{i})^{T} | — | Unit vector that points the direction of the self-driving force |

t_{d} | s | Duration for which the self-driving force is kept constant |

t_{m} | s | Maximum value of t_{d} |

Δt | s | Time step |

E_{T} | s | Time until 90% of the tasks are executed |

E_{D} | m | Total distance traveled until 90% of the tasks are executed |

W | — | Amount of remaining tasks |

D | m | Total distance traveled by all the robots |

Parameters for the improved NRIB swarm algorithm | ||

τ_{d} | s | Time constant for the change in the communication range |

d_{min} | m | Minimum value of the radius of the communication range |

d_{max} | m | Maximum value of the radius of the communication range |

γ | — | Sensitivity of the communication range to the workload |

$\alpha \u02dc$′ | 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.

*i*,

*w*

_{i}, evolves according to the following equation:

*λ*is a positive constant. Equation 2 indicates that

*w*

_{i}is given by the first-order delay of

*λz*(

**r**

_{i}) with time constant

*τ*

_{w}. Thus,

*w*

_{i}increases when robot

*i*remains in the area with high

*z*(

**r**). Note that an upper limit of

*w*

_{i}is introduced; specifically,

*w*

_{i}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 *ϵn*_{k} in unit *k*, where *ϵ* is a positive constant and *n*_{k} 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, *w*_{i} decreases as the tasks are executed, that is, *z*(**r**_{i}) approaches zero.

*i*is designed as follows:

**R**

_{ij}=

**r**

_{j}−

**r**

_{i}, $R\u02c6$

_{ij}=

**R**

_{ij}/|

**R**

_{ij}|,

*f*

_{0}is a positive constant, and

**n**

_{i}is a unit vector. The factor 1 −

*w*

_{i}is the mobility of robot

*i*. Specifically, robot

*i*slows down as

*w*

_{i}increases, and the driving force completely vanishes when

*w*

_{i}= 1. Owing to this factor, robots can continue executing tasks until they are finished. The term

*f*

_{0}

**n**

_{i}is the self-driving force, which enables random walk of the robot. The unit vector

**n**

_{i}is given by

**n**

_{i}= (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

*t*

_{d}, where

*t*

_{d}is a uniformly distributed random number in the interval (0,

*t*

_{m}] 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

*k*

^{phys}and

*μ*are positive constants and

*r*

_{0}denotes the radius of the robot.

*g*(|

**R**

_{ij}|) determines the interaction between robots. The term

*g*(|

**R**

_{ij}|)$R\u02c6$

_{ij}in Equation 3 denotes attractive and repulsive force when

*g*(|

**R**

_{ij}|) is positive and negative, respectively. We designed this function by drawing inspiration from the NRIB model [22, 23]:

*i*approaches robot

*j*, and its contribution is large when the workload of robot

*j*,

*w*

_{j}, 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*(|

**R**

_{ij}|) 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.

*α*/

*η*, $\beta \u02dc$ =

*β*/

*η*, $f\u02dc$

_{0}=

*f*

_{0}/

*η*, and $f\u02dcijphys$ =

*η*

^{−1}$fijphys$. Note that the absolute value of $r\u02d9$

_{i}is reset to

*v*

_{max}when it exceeds

*v*

_{max}.

### 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 *w*_{i} 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 **r**_{i} (Equation 6), which often varied rapidly, was solved with the fourth-order Runge-Kutta method, while that of *w*_{i} (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.

Variable . | Dimension . | Meaning . |
---|---|---|

p_{0} | — | Parameter related to the waiting time |

c | — | Parameter related to the waiting time |

n_{i} | — | Number of robots within the communication range |

τ_{bee} | s | Duration for the measurement of t_{r,i}, I_{min,i}, and I_{max,i} |

n_{col,i} | — | Number of collisions during the duration τ_{bee} |

t_{r,i} | s | Average time between subsequent collisions over the duration τ_{bee} |

$t\u02dc$_{r,i} | s | Average of t_{r,i} over the robots within S_{i} |

I_{min,i} | — | Minimum value of z(r_{i}) during the duration τ_{bee} |

I_{max,i} | — | Maximum value of z(r_{i}) during the duration τ_{bee} |

$I\u02dc$_{min,i} | — | Minimum value of I_{min,i} for the robots within S_{i} |

$I\u02dc$_{max,i} | — | Maximum value of I_{min,i} for the robots within S_{i} |

Ī | — | Normalized amount of tasks |

δ | — | Parameter introduced to avoid zero divide (Equation 10) |

Variable . | Dimension . | Meaning . |
---|---|---|

p_{0} | — | Parameter related to the waiting time |

c | — | Parameter related to the waiting time |

n_{i} | — | Number of robots within the communication range |

τ_{bee} | s | Duration for the measurement of t_{r,i}, I_{min,i}, and I_{max,i} |

n_{col,i} | — | Number of collisions during the duration τ_{bee} |

t_{r,i} | s | Average time between subsequent collisions over the duration τ_{bee} |

$t\u02dc$_{r,i} | s | Average of t_{r,i} over the robots within S_{i} |

I_{min,i} | — | Minimum value of z(r_{i}) during the duration τ_{bee} |

I_{max,i} | — | Maximum value of z(r_{i}) during the duration τ_{bee} |

$I\u02dc$_{min,i} | — | Minimum value of I_{min,i} for the robots within S_{i} |

$I\u02dc$_{max,i} | — | Maximum value of I_{min,i} for the robots within S_{i} |

Ī | — | Normalized amount of tasks |

δ | — | Parameter introduced to avoid zero divide (Equation 10) |

Parameter . | Value . | ||
---|---|---|---|

Condition 1 . | Condition 2 . | Condition 3 . | |

Nondimensional unit length | 1 m | 2 m | 3 m |

Nondimensional unit time | 1 s | 2 s | 3 s |

Parameter . | Value . | ||
---|---|---|---|

Condition 1 . | Condition 2 . | Condition 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 *t*_{maxstep}.

First, simulations were performed with various *d*_{0} 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 *d*_{0} 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 *d*_{0} = 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 *d*_{0} = 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 *d*_{0} = 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)).

*E*

_{T}and

*E*

_{D}, respectively. Namely,

*E*

_{T}and

*E*

_{D}are defined as

**v**

_{i}(

*k*) is the velocity vector of robot

*i*at the

*k*th time step and

*T*is the time step when 90% of the tasks in the field are executed, that is, ∫

_{F}

*z*(

**r**)dr decreased by 90% from the initial condition, with ∫

_{F}being the spatial integral over the field. Note that

*T*is set to

*t*

_{maxstep}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 *E*_{T} and *E*_{D} are the smallest when *d*_{0} is 14 and 6 m, respectively. Thus, the robots can execute tasks fast with a short travel distance when *d*_{0} is around 10 m.

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

*t*_{bee,i}, which is set to become longer as the amount of tasks*z*(**r**_{i}) 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 *t*_{bee,i} is longer as the amount of tasks *z*(**r**_{i}) 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.

*t*

_{bee,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*(

**r**

_{i}) 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

*t*

_{bee,i}is designed as

*p*

_{0}and

*c*are positive constants, $t\u02dc$

_{r,i}denotes the average robot-to-robot encounter time interval, and $I\xaf$

_{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 *E*_{T} and *E*_{D} values when *p*_{0} is varied. A representative trial when *p*_{0} = 5.0 is shown in Movie 4 (see the online supplementary material). For all *p*_{0} values, *E*_{T} and *E*_{D} are larger than those for the proposed control scheme. According to the statistical analysis, the *E*_{T} and *E*_{D} values when *d*_{0} = 10 m in the proposed control scheme are significantly smaller than those when *p*_{0} = 5 in the BEECLUST algorithm (Welch test, *p* < 0.01).

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 *d*_{0} is varied are shown in Figure 4. Both *E*_{T} and *E*_{D} decrease as *d*_{0} increases, and do not vary significantly for large *d*_{0}. Thus, the optimal value in Condition 1 (*d*_{0} ≃ 10 m, Figure 2) is not optimal in conditions 2 and 3. The optimal value of *d*_{0} 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.

## 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 *d*_{i} 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.

*d*

_{i}is designed as

*γ*and

*τ*

_{r}are positive constants. The radius of the communication range

*d*

_{i}approaches

*d*

_{max}and

*d*

_{min}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.

*d*

_{i}is large, robot

*i*is required to approach distant robots with high workload to execute tasks cooperatively. Thus, the attraction term ($\alpha \u02dc$

*w*

_{i}|

**R**

_{ij}|

^{−1}in Equation 6) should increase as

*d*

_{i}increases. Hence, Equation 6 is modified as

### 5.1 Simulation Results

Simulations were performed using the improved NRIB swarm algorithm. The time evolution of *d*_{i} (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 *E*_{T} and *E*_{D} 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, *E*_{T} and *E*_{D} for the improved NRIB swarm algorithm are significantly smaller than those for the NRIB swarm algorithm with *d*_{0} = 10 m (Welch test or *t*-test, *p* < 0.01), except for *E*_{D} for condition 1, wherein no significant difference was found. This was achieved by auto-tuning of the communication range.

Parameter . | Value . | |||||
---|---|---|---|---|---|---|

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 |

t_{maxstep} | 2000000 | 2000000 | 2000000 | 2000000 | 2000000 | 2000000 |

r_{0} | 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 |

λ | 5 | 5 | 5 | 5 | 5 | 5 |

ϵ | 0.1 | 0.1 s^{−1} | 0.2 | 0.1 s^{−1} | 0.3 | 0.1 s^{−1} |

$\alpha \u02dc$ | 2.000 | 2 m^{2} s^{−1} | 1.000 | 2 m^{2} s^{−1} | 0.667 | 2 m^{2} s^{−1} |

$\beta \u02dc$ | 2.000 | 2 m^{3} s^{−1} | 0.500 | 2 m^{3} s^{−1} | 0.222 | 2 m^{3} s^{−1} |

$f\u02dc$_{0} | 0.5 | 0.5 m · s^{−1} | 0.5 | 0.5 m · s^{−1} | 0.5 | 0.5 m · s^{−1} |

$k\u02dc$_{phys}(≡ k_{phys}/η) | 1000 | 1000 m^{−4} s^{−1} | 32000 | 1000 m^{−4} s^{−1} | 243000 | 1000 m^{−4} s^{−1} |

μ | 5 | 5 | 5 | 5 | 5 | 5 |

v_{max} | 1 | 1 m · s^{−1} | 1 | 1 m · s^{−1} | 1 | 1 m · s^{−1} |

t_{m} | 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 |

Parameter . | Value . | |||||
---|---|---|---|---|---|---|

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 |

t_{maxstep} | 2000000 | 2000000 | 2000000 | 2000000 | 2000000 | 2000000 |

r_{0} | 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 |

λ | 5 | 5 | 5 | 5 | 5 | 5 |

ϵ | 0.1 | 0.1 s^{−1} | 0.2 | 0.1 s^{−1} | 0.3 | 0.1 s^{−1} |

$\alpha \u02dc$ | 2.000 | 2 m^{2} s^{−1} | 1.000 | 2 m^{2} s^{−1} | 0.667 | 2 m^{2} s^{−1} |

$\beta \u02dc$ | 2.000 | 2 m^{3} s^{−1} | 0.500 | 2 m^{3} s^{−1} | 0.222 | 2 m^{3} s^{−1} |

$f\u02dc$_{0} | 0.5 | 0.5 m · s^{−1} | 0.5 | 0.5 m · s^{−1} | 0.5 | 0.5 m · s^{−1} |

$k\u02dc$_{phys}(≡ k_{phys}/η) | 1000 | 1000 m^{−4} s^{−1} | 32000 | 1000 m^{−4} s^{−1} | 243000 | 1000 m^{−4} s^{−1} |

μ | 5 | 5 | 5 | 5 | 5 | 5 |

v_{max} | 1 | 1 m · s^{−1} | 1 | 1 m · s^{−1} | 1 | 1 m · s^{−1} |

t_{m} | 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 *t*_{maxstep} 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 *d*_{0} was set to 10 m in the NRIB swarm algorithm, and *p*_{0} was set to 5.0 in the BEECLUST algorithm. The performance was evaluated by *E*_{T} and *E*_{D} in cases (ii) and (iii); however, *E*_{T} could not be used in case (i) because of the appearance of tasks. Hence, the performance in case (i) was evaluated by *E*_{D} and *E*_{R}, the latter of which is defined as the ratio of ∫_{F}*z*(**r**)d**r** at *t* = *t*_{maxstep} to that at *t* = 0.

The results are shown in Figure 5. In case (i), *E*_{R} and *E*_{D} for the BEECLUST algorithm are significantly larger than for the other algorithms (Welch test, *p* < 0.01). *E*_{R} for the improved NRIB swarm algorithm is significantly smaller than that for the NRIB swarm algorithm (Welch test, *p* < 0.01). However, *E*_{D} 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), *E*_{T} and *E*_{D} 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).

## 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 *d*_{0} is an important parameter. For small *d*_{0}, the robots move almost independently and need to travel for a long time to find tasks, while for large *d*_{0} they move long distances to aggregate (Figure 1). Thus, *d*_{0} 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*(**r**_{i}) 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

### 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.

Basically, each robot moves with the velocity vector v_{max}**n**_{i}, where **n**_{i} = (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).

*t*

_{r,i},

*I*

_{min,i}, and

*I*

_{max,i}defined below and shares them with other robots within its communication range.

*t*_{r,i}: Average time between subsequent collisions over the duration*τ*_{bee}in the recent past. It is defined by*t*_{r,i}=*τ*_{bee}/*n*_{col,i}, where*n*_{col,i}denotes the number of collisions during the duration*τ*_{bee}.*I*_{min,i}: Minimum value of*z*(**r**_{i}) during the duration*τ*_{bee}in the recent past.*I*_{max,i}: Maximum value of*z*(**r**_{i}) during the duration*τ*_{bee}in the recent past.

_{r,i}, $I\u02dc$

_{min,i}, and $I\u02dc$

_{max,i}given by

*n*

_{i}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\u02dc$

_{i}is calculated as

_{r,i}and $I\xaf$

_{i}, the waiting time

*t*

_{bee,i}is given by Equation 9 in the main text.

When the robot collides with another robot, it stops motion for the duration *t*_{bee} and then moves with the velocity vector *v*_{max}**n**_{i} 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.