Estimating the Expected Time to Enter and Leave a Common Target Area in Robotic Swarms

1. Introduction

1.1. Motivation

Robotic swarms are systems formed by often simple robots with distributed processing and local sensing [1]. As it happens in insect swarms, no centralised coordination is performed for controlling each individual. Thus, distributed control and local sensing are assumed in robotic swarms [2,3].

As a way to reduce costs, robotic swarms use many simple and cheap robots instead of relying on a few complex and expensive ones [4]. Creating complex group behaviours based on simple rules is the goal of some robotic swarm research [1,5]. Counter-intuitively, such complex global behaviours may emerge from well-planned local rules. [6,7]. Consequently, swarm designers have to infer from local rules a global result needed to solve their problems. This inference may arise from several simulations of the local rules. However, when dealing with many robots, simulations could be time-consuming, and estimations of the resulting global outcome would be useful for swarm designers.

Having that in mind, this paper focuses on estimating a global outcome, the completion time of a task performed by some number of robots in a swarm, based on the robot’s local specification in a problem that may occur when numerous robots of that swarm have the same area: the common target problem [8]. This problem may occur, for instance, when robots are deployed to collect a leaked toxic substance. These robots need to coordinate themselves to reach that region as fast as possible while not interfering with the collection of the other robots. If the robots do not execute an efficient traffic control algorithm, their mutual interference can worsen.

Hence, solutions to the common target problem can also be applied in waste collection [9,10] in secluded areas such as the ocean, disaster response [11], construction [12], searching of cancer cells by microscopic robots [13,14,15,16], rescue in regions harmful to humans [17] and warehouse management [18,19]. In this problem, the number of robots inevitably may cause congestion as they are forced to share the same area. In addition, this number can impact the execution time, which affects other measures (e.g. heat generation and energy consumption [20]).

Thus, when devising different solutions for this problem, the typical practice is to compare the task completion time with the number of robots and environmental or algorithmic parameters (e.g., working area and target area sizes) [8,21,22,23]. For the common target problem, artificial potential fields are usually applied to such solutions in the same manner as in robotic swarms due to their distributed and local characteristic for path planning [24,25,26]. Additionally, other collision methods may not work well for this problem. For instance, Marcolino et al. [8] performed tests with the distributed collision avoidance method ORCA, which does not work for the common target problem because the robots enter a circular deadlock, not reaching the target area. The recent GBP Planning [27] requires the additional cost of using communication between robots, while the algorithms mentioned in this work do not need it.

That being said, when the local individual specification based on potential fields is modified, global swarm behaviours, such as the task completion time, will have specific outcomes accordingly. To compare distinct potential fields that work better for various numbers of robots, experiments with different values must be simulated to form a plot of the function of the task completion time versus the number of robots or other parameters. Nevertheless, plotting it for many robots through simulation is time-consuming. A swarm designer could use estimations instead of simulations to choose between available solutions.

An alternative approach would be using regression over the experimental data. However, it does not explain the time in relation to the number of robots, dynamical variables such as average speed, the average distance between the robots, and environmental or algorithmic variables such as the radius of the target area, radius of the task working area and how many directions the robots must go after they reach the target area.

Additionally, there are many situations where different swarm teams may be present in the same environment. For instance, such mixed teams (MT) can be found in the multi-agent literature [28] in problems like Pursuit-Evasion Games [29,30] and Multi-Agent Pickup and Delivery [31]. This work considers mixed teams where each group does not know the algorithm used by the other groups. This could happen, for example, in a situation where a robotic swarm runs some coordination algorithm for accessing an area for rescuing people in an accident, and different institutions want to help. Assume that they deployed their robots and, for a fast response, did not meet to settle the swarm congestion control algorithm for quickly accessing the common area. Then, the first swarm to arrive will execute a congestion control algorithm, but the next does not know which algorithm is in execution.

This is an example of ad hoc teamwork [32] and usually agents use learning about other agents to handle it [33,34,35]. However, the above situation requires a fast response, but learning needs long iterations of training, and errors may be introduced during the learning process. Thus, it is considered here that each group of robots in MT executes a fixed algorithm without spending time on learning the algorithm of others. If these algorithms happen to be the same across all groups there will be no impact in the swarm by mixing them. However, if they are different, they may affect each other negatively. Therefore, the influence of different algorithms is also analysed here. In fact, this work shows that the expected task completion time estimations are also applicable to estimate the expected task completion time of mixed teams.

1.2. Objective

The objective of this paper is to introduce estimations of the expected task completion time according to the number of robots and environmental and algorithmic parameters in algorithms for the common target problem using potential fields. As these estimations are difficult to be expressed directly from the individual robot controller, they are based on analysing the observations made by simulating the control swarm algorithm. The equations of these estimations are mainly derived from geometric characteristics of the usual behaviour of the robots by following the control algorithm or obtained from well-known equations cited in the explanation of the estimation. No paper has presented such estimations for that problem before. Thus, this paper is a fundamental step to initiate a discussion about how the global system behaviour can be inferred from individual potential field-based local controllers in Swarm Robotics.

For doing so, extensive experiments were performed in the realistic Stage simulator [36] due to its easiness of running in batch and parallel, low memory cost, fast execution for numerous robots and recent usage and review [37,38,39], assessing the completion time from hundreds of robots as quick as possible to generate graphs and analyse different parameters. The proposed estimations are compared with data from experiments for a swarm with no coordination [21], state-of-the-art algorithms Single Queue Former and Touch and Run Vector Fields [23] and two groups in MT. Also, these algorithms are modelled in equation form to highlight the difficulty of directly deducing the task completion time from them. Experiments with them were made from 20 to 300 robots and different circular target area radii with holonomic and non-holonomic robots, to show the adaptability of the estimations to divergent types of robot movement. For MT, all combinations of those algorithms were tested to show the suitability of each algorithm for different ratios of the number of robots in a group to the total number of robots. The root mean squared error normalised by the y-axis range (NRMSE) values of the estimations are between 0.011 and 0.028 for a swarm without different groups and between 0.008 and 0.036 for two groups in MT.

1.3. Paper Organisation

This paper is organised as follows. The next section discusses the related work, including studies of local-to-global in science, the many-body problem and its application in Swarm Robotics and the robotic swarm common target problem and the algorithms to solve it.

Section 3 presents the task completion time estimations and how they were obtained for a swarm using only attractive and repulsive force fields, state-of-the-art algorithms and Mixed Teams. In that section, each algorithm is stated in its equation form and exemplified with key screenshots of some executions. These examples are samples of the observations used for describing how these estimations are approximated. Thus, following each estimation, an explanation shows the derivations of these approximations developed from those observations made by simulating the algorithm.

Section 4 describes the experiments performed to corroborate the estimations derived in this work and discusses their results. These experiments aim to compare the estimations with the data obtained from simulations.

Section 5 discusses the limitations, challenges, improvements over this work and how the methodology can be applied to similar problems in Swarm Robotics.

Finally, Section 6 summarises the results and concludes this work.

2. Related Work

The effect of the task completion time of a robotic swarm based on the dynamic of one individual is an example of an analysis from local to global. Similar works can be found in other study areas. In special, the many-body problem in Physics has elegant and strong mathematical solutions. Some works for traffic problems or Swarm Robotics have applied solutions for the many-body problem. However, the solutions presented are not adequate for the common target problem. The algorithms developed for the common target problem mostly use potential fields as a local path planner to emerge global congestion control. These works will now be discussed in detail.

Studies of Local-to-Global: Biology, Chemistry and Physics share the study of individual parts of different systems to explain their global behaviour. The manifestation of emergence [40,41,42,43] is a good example of such studies. In it, the microscopic characteristics enable the arising of unexpected macroscopical behaviour. As the understanding of these relations grows, advances in other areas like engineering and robotics are made. Related to that, the interdisciplinary Synergetics [44] explain self-organised phenomena in systems with non-equilibrium such as cloud dissipation, lasers and patterns in slime mould or chemical reactions.

The Many-Body Problem and Solutions: In Physics, problems about the interaction of microscopic particles is generally mentioned as the many-body problem [45,46]. For instance, in statistical mechanics, the Maxwell-Boltzmann distribution function [47,48] states the speed probability distribution of ideal gas particles. It is calculated by analysing the movement of the particles under the influence of gravity and the probability frequency of the energy from the distance to the ground. Then, the velocities are averaged over all directions and that probability frequency to yield the desired distribution. There are different approaches to this problem. The perturbation theory [49,50] starts with a simpler version of the problem. Then, an expression, often called perturbative expansion, using power series of successive corrections to the simpler problem is developed to solve the original problem. The coupled cluster technique is used in computational chemistry [51,52] for solving this problem numerically. Similarly, the density-functional theory [53,54] is applied in computational modelling to understand the electronic and molecular organisation from atoms to solids by using functionals of electron density that depend on the space. The mean-field theory [55] substitutes the high number of interactions between the particles with a few averaged ones, reducing the number of calculations but also the exactness of the solution. There are interesting applications of this method in multi-agent research such as Chen et al. [56], Lasry and Lions [57]. The variational method [58,59] is a solution similar to that used in quantum mechanics for approximating the desired energy state by selecting wave functions whose parameters minimise the expected value of the energy. The neural network quantum states [60] approximate the solution using machine learning to minimise the variational energy.

Similar to Many-Body in Swarm Robotics: Additionally, the literature on many-body problem is mostly Physics-oriented. An interesting approach to the intersection of Physics and Swarm Robotics is described by Hamann [61], where he applied Brownian motion, Langevin equation and Fokker-Plank equation for estimating global behaviour from individual robot specification by random movement. Another approach to the above solutions for the many-body oriented to real-world application is the Herman-Prigogine kinetic equations for vehicular traffic [62,63,64], considering the drivers with intentions but random actions. However, as these approaches give results assuming that every action of every driver or particle is random as time passes, their equations cannot be applied directly in this work because only the initial position is random, but all the movement is deterministic. The survey of Elamvazhuthi and Berman [65] also presents another study on Swarm Robotics and Physics aligned with the problem in this paper, but the references provided by Elamvazhuthi and Berman [65] do not directly answer the research question shown here. Instead of assuming randomness in the movement of the individual robots and showing purely theoretical deductions, the estimations calculated here use experimental observation from simulations performed in Stage [36]. Additionally, these estimations are a fundamental first step to start a discussion on theoretical aspects on the relation of time and the number of robots on the common target problem.

The Robotic Swarm Common Target Area Problem and Algorithms: The common target problem is usually handled by using potential fields [8,21,23]. In fact, the algorithms discussed in Section 3 use them in the role of an individual path planner to accomplish global congestion control. Any of these algorithms modify the attractive and the repulsive force field to minimise the time to access the common target area or the task completion time, which, in addition to the time to access it, includes the time for leaving it. This may bring to light the question about whether it is possible to find the exact attractive potential field that minimises this time by a method similar to the Lagrangian in mechanics [66], which elegantly solves curves with minimum energy requirement using stationary action or with variability as in the many-body approaches discussed above. For finding such a minimal potential field, the function of the task completion time per number of robots has to be obtained beforehand. As far as the authors know, no work presents how to obtain this function from the motion equations of a robot, which is a challenge in Swarm Robotics. Thus, experimental approximations are presented in this paper by observing the swarm, making this work a fundamental first step for studying the common target problem as a local-to-global analysis to serve as a basis for comparing new approaches for a theory capable of finding such minimal potential field.

3. Expected Task Completion Time Estimation Analysis

In this section, the estimation analysis is presented by describing the task scenario of the experiments, the individual robot controller, the challenge in the development of a global analysis from local controllers and the common settings of the algorithms for which the expected task completion time estimations are analysed for a number of robots in the swarm.

Task Scenario: Consider the task scenario where a system with N robots in a plane must reach a common target located at $\mathbf{G}=(G_x,G_y)$. Each robot knows its position and the positions of the targets, and it can sense nearby objects (for instance, by using laser scanners covering 360 degrees). Upon reaching the target, the robots move to different destinations, which may not be similar. Suppose that the target area is a circle of radius s, and a robot reaches the target if its centre of mass is at a distance below or equal to the radius s from the centre of the target. In addition, there is no minimum amount of time to stay at the target.

In this scenario, the coordination algorithms, if applied, are only employed inside a working area, that is, a circle of radius $D>s$ around the target. Although the swarm density is not considered as an explicit parameter in this work, the number of robots N and the radius D of the working area yield the swarm density $N/\left(\pi D^2\right)$. If the robots have two or more targets, the targets are apart by at least $2D$. The initial position of the robots has a distance from the target centre randomly chosen from a uniform distribution in $[D,D+E$ m for a fixed E. The angle between the initial position of a robot and the x-axis is selected from a uniform distribution in $[0,2\pi$ rad. To avoid robots starting in the same position and crashes between robots in the beginning of the simulation, if the distance between a robot and any other is less than 1 m, its initial position is selected again. This arrangement inside an annulus is only considered at the start of an experiment and before the initial target.

Once the robots have reached the first target area, they will proceed to the next one at a long distance, either to the left or right of the shared target. In other words, the next target is chosen randomly between positions $(-999999,G_y)$ and $(999999,G_y)$. Thus, approximately half of the robots will go to the left-hand side and the rest to the right, according to a uniform probability. This paper focuses on the congestion around the first target area, which is shared by all the robots in the swarm. After they leave the working area around the shared target area, they head towards the next target area by using only repulsive and attractive forces (that is, they turn off the coordination algorithm). In addition, the robots measure $0.44\times 0.44\times 0.44$ ${\mathrm{m}}^3$ in the experiments. Although some settings are precise (e.g., robot size and next target positions), the equations derived in this paper are based on general behaviours and are applicable with different parameters.

The global metric analysed in this paper from the local individual robot controller is defined below.

Definition 1.

The task completion time is obtained from the last robot to leave the working area around the shared target area, that is, the maximum time to enter the shared target area plus the time to leave the working area for every robot in the swarm.

Once a robot leaves the working area after it reached the first target area, it is no longer considered in the computation of any metric presented in this paper. As the initial position of the robots is random in an experiment and influences the variation of the task completion time, it is considered the expected task completion time from various experiments, which is defined below.

Definition 2.

Let $t_{Alg}$ be the task completion time of an algorithm $Alg$ using N robots for an experiment with the initial position of the N robots fixed. Then, $\mathbb{E}\left[t_{Alg}\right]$ stands for the expected task completion time over all executions for N robots.

Individual Robot Controller: Consider a robot $i\in \{1,\dots ,N\}$ at position ${\mathbf{r}}_i\left(t\right)=(r_{ix}\left(t\right),r_{iy}\left(t\right))$ at time t. In this work, each robot uses the same movement controller, chosen to physically describe the change in the position of the robot given by any algorithm that has an input force resultant from a potential force field $U\left({\mathbf{r}}_i\left(t\right)\right)=U_A\left({\mathbf{r}}_i\left(t\right)\right)+U_R\left({\mathbf{r}}_i\left(t\right)\right)$, for an attractive potential field $U_A\left({\mathbf{r}}_i\left(t\right)\right)$ and a repulsive potential field $U_R\left({\mathbf{r}}_i\left(t\right)\right)$. The input force has to be equal among the robots for all algorithms to perform a fair comparison among different approaches. Also, the controller needs to adapt the acceleration when the robot’s velocity changes in situations such as the sudden advent of another robot or the transition to a new target. To satisfy these requirements, the controller is based on the previous works on the common target problem [8,21,23] but rewritten in continuous form: $\dot{{\mathbf{r}}_i}\left(t\right)={\displaystyle \frac{d{\mathbf{r}}_i\left(t\right)}{dt}},\ddot{{\mathbf{r}}_i}\left(t\right)={\displaystyle \frac{d\dot{{\mathbf{r}}_i}\left(t\right)}{dt}}$, and

$$\ddot{{\mathbf{r}}_i}\left(t\right)=K_{res}{\displaystyle \frac{{\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)+{\sum }_{\begin{array}{c}j=1\\ i\ne j\end{array}}^N{\mathbf{F}}_{\mathbf{R}}({\mathbf{r}}_i\left(t\right),{\mathbf{r}}_j\left(t\right))}{∥{\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)+{\sum }_{\begin{array}{c}j=1\\ i\ne j\end{array}}^N{\mathbf{F}}_{\mathbf{R}}({\mathbf{r}}_i\left(t\right),{\mathbf{r}}_j\left(t\right))∥}}-K_{dp}\dot{{\mathbf{r}}_i}\left(t\right),$$

(1)

where ${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)=-\nabla U_A\left({\mathbf{r}}_i\left(t\right)\right)$ is a force to attract the robot toward a target, and ${\mathbf{F}}_{\mathbf{R}}\left({\mathbf{r}}_i\left(t\right)\right)=-\nabla U_R\left({\mathbf{r}}_i\left(t\right)\right)$ is the repulsive force applied to avoid bumping other robots. Without the assumption that the robot knows the target position, ${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)$ cannot be computed correctly. ${\mathbf{F}}_{\mathbf{A}}$ depends on the algorithm used by each robot and has a fixed magnitude $K_a$. The sections below show different choices for ${\mathbf{F}}_{\mathbf{A}}$. Each robot adds the attractive force with the repulsive forces and then constrain the result to a fixed vector modulus $K_{res}$ to homogenise the acceleration among the individuals in the swarm. The constant $K_{dp}>0$ is employed to control alterations in velocity. It prevents the robots slipping due to acceleration changes.

Also, the following repulsive force is employed because it acts as a strong barrier as the robot approaches an obstacle [67]:

$${\mathbf{F}}_{\mathbf{R}}({\mathbf{r}}_i\left(t\right),{\mathbf{r}}_j\left(t\right))=\left\{\begin{array}{cc}{K}_r\left({\displaystyle \frac{1}{\parallel {\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j\left(t\right)\parallel }}-{\displaystyle \frac{1}{I_d}}\right){\displaystyle \frac{{\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j\left(t\right)}{\parallel {\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j{\left(t\right)\parallel }^3}},\hfill & \mathrm{if}\parallel {\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j\left(t\right)\parallel <I_d\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

where $K_r$ is a fixed multiplicative constant for the repulsive force field, and $I_d$ is the default influence radius, that is, the maximum distance from its mass centre a robot considers anything sensed as an obstacle to avoid. $I_d$ can be set to the maximum range of the robot’s local sensor. Thus, a robot deprived of local object sensing cannot compute ${\mathbf{F}}_{\mathbf{R}}({\mathbf{r}}_i\left(t\right),{\mathbf{r}}_j\left(t\right))$.

Challenge in the Development of a Global Analysis from Local Controllers: With that ${\mathbf{F}}_{\mathbf{R}}$, Equation 1 depends on the neighbourhood of the robot, which, although limited, is dynamic. Even if there is a fixed maximum number of neighbours, a global analysis of ${\mathbf{r}}_i\left(t\right)$ would depend on all the other $N-1$ robots for each instant t. Thus, obtaining a closed-form expression for the task completion time in terms of N directly from Equation 1 is difficult. In addition, observing the behaviour of an algorithm for hundreds of robots in the real world is expensive, but it is faster and more affordable by simulation.

In other words, if the individual robot control equations were simpler than Equation 1, the closed-form expression of the expected task completion time in terms of N would be obtained by the following steps: (i) integrate Equation 1 twice to obtain ${\mathbf{r}}_i\left(t\right)$ for $i\in 1,\dots ,N$, (ii) invert it to retrieve t for each robot i, (iii) obtain the maximum t for every robot i and (iv) calculate the expected value of this maximum when the initial positions of the robots are randomly distributed inside the annulus of radii D and $D+E$ and centre at $\mathbf{G}$. As this closed-form expression cannot be derived from Equation 1 in such manner, a theoretical macroscopic analysis inspired by simulations is performed to attain an approximation of the expected task completion time $E\left[t\right]$. For doing this, besides the target area radius s and the working radius D, the following variables are given: the mean distance between the centre of mass of a robot and the others $\overline{d}$ inside the influence radius $I_d$ for all robots and the mean linear speed of all robots $\overline{v}$. A summary of the variables most used henceforward is in Table 1.

Moreover, the approximations stated in this section may use these variables as input and are mainly derived from geometric observations or obtained from well-known equations. These approximations are explained after their statements, based on a macroscopic analysis observed from the simulations. As these approximations are constructed from experimentation, these explanations cannot be considered formal proofs, although they are composed of proven mathematical expressions.

Furthermore, in the expected task completion time expressions shown below, constants are defined for each algorithm to subsume the overall behaviour of the robots, their interactions and the congestion inspired by the experiments. Such abstraction is analogous to the friction constant in Physics, which summarises the microscopical effect of surfaces. In this field of study, the friction constant for commonly used materials was initially calculated from experiments with them. Similarly, the defined constants are fitted by the experimental data in Section 4.

Common Settings of the Presented Algorithms: In this section, four algorithms are analysed about the expected task completion time. In the examples in the figures that follow and in the performed experiments, it was used the values in Table 2 and the robots are holonomic. Note, however, that the theoretical results also apply for arbitrary values and holonomic or non-holonomic robots, as shown in Section 4.

In the following sections, the state-of-the-art algorithms in [23] are restated in equation form and based on the local motion controllers for showing the difficulty of solving the task completion time from the individual local control equations. This presentation begins with the no coordination (NC) [21] (that is, without any congestion control algorithm, only attractive potential field to the target area and repulsive potential field for avoiding other robots) as a simple example to show that, even from this reduced model, it is hard to infer the expected task completion time per number of robots.

3.1. No Coordination Algorithm

Algorithm Description: When the robots are not coordinated, they only follow the target region, avoiding the others by repulsive force. The equation for the ${\mathbf{F}}_{\mathbf{A}}$ in Equation 1 is the attractive force to the target centre $\mathbf{G}$ given by

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)={\displaystyle \frac{K_a(\mathbf{G}-{\mathbf{r}}_i\left(t\right))}{\parallel \mathbf{G}-{\mathbf{r}}_i\left(t\right)\parallel }}.$$

(2)

Although this is the simplest attractive force of all algorithms presented here, Equation 1 also considers the repulsive force of all robots in the environment. Thus, as the number of robots grows, the analytical solution for the differential equations is still complicated to infer.

Samples of the Observations: As an illustration, Figure 1 and Figure 2 show an experiment of NC. Red robots are going to the target region, yellow ones have arrived at it and are trying to leave it, and black robots have exited from the working area and are going to the next target. The robots start in random positions outside the working area (Figure 1(a)). The first robot to reach the working area takes an expected time inversely proportional to the number of individuals in the swarm (marked with a circle and an arrow in Figure 1(b)). Each robot goes toward the target area until the first reaches it . Eventually, the first robots to reach the common target area will be stuck in it and will try to leave while they are slowed down due to the congestion caused by the other robots that did not arrive at the target yet (Figure 1(c)). The robots trying to leave the target area will slowly push the other robots until they move out of the cluttered region with a shape similar to a circle. If a robot is heading in the same direction as the first to overcome this area, it may follow the space created by the other robots avoiding the first. This may cause more robots to do the same, and a queue of robots is formed (Figure 1(d)). The length of this queue and the number of robots in it are random. New queues may appear while the number of robots in the cluttered area diminishes. When no robots are trying to arrive at the target area, the last queues of robots will finally reach the outside of the working area (Figure 2(e)) until the last robot leaves this area (with a circle and an arrow in Figure 2(f)).

Figure 3 and Figure 4 illustrate an example with a larger target area radius. The robots have more space to access the target. In this example, the capacity of the target area is greater than the number of robots, so they can use the free space to travel to the next targets in both sides (Figure 3(c)). Due to this free space, the number of robots going to the target (in red) is not enough to force the robots leaving it (in yellow) to go to the target centre (Figure 3(d)). Thus, the time until the cluttering near the target area disappears is lower than in the previous example, and how the robots leave the target area is also different (Figure 4(e)) because of the free space in the target area and the less resistance caused by the robots going to the target area.

By analysing such examples, an approximation of the expected task completion time for the NC algorithm can be obtained. The following estimation presents the result. After its statement, an explanation shows the derivations of that approximation based on the observations made by simulating the algorithm. Because the following estimation is constructed from experimentation, the explanation below cannot be considered a formal proof.

Estimation 3.

The estimated expected time for N robots starting at a random position with a distance from the target centre in $[D,D+E$ to arrive at the common target area and leave the working area without coordination is

$$\mathbb{E}\left[t_{NC}\right]\approx {\displaystyle \frac{E}{\overline{v}(N+1)}}+{\displaystyle \frac{2D-s}{\overline{v}}}+wai{t}_1\left(N\right).$$

(3)

for

$$wai{t}_1\left(N\right)=\left\{\begin{array}{cc}{C}_{NC1}N,\hfill & \mathrm{if}N\le {\displaystyle \frac{2{\left(s+\overline{d}/2\right)}^2}{{\left(\overline{d}/2\right)}^2}}\hfill \\ C_{NC2}\left({\displaystyle \frac{N^{1.5}}{1.5\sqrt{\pi }}}-N\right),\hfill & \mathit{otherwise}\hfill \end{array}\right.$$

and constants $C_{NC1}$ and $C_{NC2}$.

Explanation. 

The equation for estimating the expected time for NC has five parts: the expected time for (a) the first robot arriving at the working area from its starting position, (b) the first robot that entered the working area reaching the target area (e.g., the time at Figure 1(b)), (c) the first robots filling the target area (e.g., the time from Figure 1(b) to 1(c)), (d) the cluttered area disappearing by the queues of robots leaving the target (e.g., the time of Figures 1(c)–2(e) and 3(c)-4(e)) and (e) the last robot that reached the target leaving the working area (e.g., the time from Figure 2(e) to 2(f)).

Part (a) is obtained from the expected value of the minimum starting distance to the working area for the N robots. Due to the task scenario description, this value uniformly ranges in $[0,E$. Consequently, the expected minimum distance is $E/(N+1)$, obtained from the probability distribution function of the minimum value [68, p. 229]. Then, the expected time for the first robot to arrive at the working area is $E/\left(\overline{v}(N+1)\right)$.

Parts (b) and (e) share the same expression, obtained from the expected time for a robot on average speed going from the end of the working area to the target area and vice-versa: $(D-s)/\overline{v}$. Part (c) is given by the expected time for a robot going from the border of the target area to the target centre due to the congestion caused by the other robots going to the target area: $s/\overline{v}$.

For calculating the time (d), there are two cases. Let the capacity of the target area be the maximum number of robots that fit inside the target area. When the capacity of the target area is greater than or equal to the number of robots outside the target area, the robots which arrive at the target area have less resistance to leave through the space between the robots which did not yet. The capacity of the target area is approximately the number of circles of radius $\overline{d}/2$ – half the mean distance between the robots – inside the target area with radius s increased by $\overline{d}/2$, that is, ${\displaystyle \frac{{\left(s+\overline{d}/2\right)}^2}{{\left(\overline{d}/2\right)}^2}}$. Thus, the number of robots exceeding the capacity of the target area is less than or equal to its capacity if $N-{\displaystyle \frac{{\left(s+\overline{d}/2\right)}^2}{{\left(\overline{d}/2\right)}^2}}\le {\displaystyle \frac{{\left(s+\overline{d}/2\right)}^2}{{\left(\overline{d}/2\right)}^2}}\iff N\le {\displaystyle \frac{2{\left(s+\overline{d}/2\right)}^2}{{\left(\overline{d}/2\right)}^2}}$. From experiments, the waiting time can be approximated by $C_{NC1}N$ for a constant $C_{NC1}$ that abstracts the influence of other factors in this time, assuming that it does not depend on time (e.g., any scaling in the speed or the average number of robots per queue, and the influence in the overall movement by the type of the robot – holonomic or non-holonomic).

If the maximum number of robots that fit inside the target area is less than the number of robots outside the target area, the resistance for leaving the target area is higher than before, and part (d) is calculated differently. Due to the erratic formation of the queues of the robots leaving the target area, the time is approximated by assuming a fixed mean distance $\overline{d}$ between the robots and mean speed $\overline{v}$ moving from the cluttered area and taking the same time to move as the erratic queues. As seen in the example from Figure 5, the robots outside the target area go towards the centre, while those leaving it form queues following the other robots for each direction to overcome the congestion. The red robots going to the target tend to concentrate inside a circle containing the cluttered area. On the other hand, the yellow robots leave that area one by one for each queue, diminishing the number of robots in the clutter.

Hence, the cluttered area is approximated by a circle whose radius $r\left(t\right)$ vanishes as the robots run through these queues. Let $N\left(t\right)$ be the number of robots in that circle at time t. Figure 5 illustrates this circle with robots. The cross mark in the middle indicates the target area centre. Suppose that a robot next to that position has to leave that circle. Then, it must run a distance of at most $r\left(t\right)$.

This radius $r\left(t\right)$ can be estimated from the area occupied by $N\left(t\right)$ robots. For this estimation, consider each robot occupying a squared area of ${\overline{d}}^2$ ${\mathrm{m}}^2$. The area of the circle containing the robots is approximately the area occupied by $N\left(t\right)$ squares with area ${\overline{d}}^2$, that is,

$$\pi r{\left(t\right)}^2\approx {\overline{d}}^{2}N\left(t\right)\Rightarrow r{\left(t\right)}^2\approx {\displaystyle \frac{{\overline{d}}^{2}N\left(t\right)}{\pi }}\Rightarrow r\left(t\right)\approx \sqrt{{\displaystyle \frac{{\overline{d}}^{2}N\left(t\right)}{\pi }}}\Rightarrow r\left(t\right)\approx \overline{d}{\displaystyle \frac{\sqrt{N\left(t\right)}}{\sqrt{\pi }}}.$$

A robot next to the target centre has to move through a distance equivalent to the number of robots fitting $r\left(t\right)$, excluding that one. Thus, the number of robots to move through is given by $r\left(t\right)/\overline{d}-1$. As one robot is decreased proportionally to the number of robots that it would have to move through and the number of queues, the rate of decrease is given by

$${\displaystyle \frac{dN\left(t\right)}{dt}}={\displaystyle \frac{-1}{C_{NC2}\left(\frac{r\left(t\right)}{\overline{d}}-1\right)}}=-{\displaystyle \frac{1}{C_{NC2}\left(\frac{r\left(t\right)}{\overline{d}}-1\right)}}\approx -{\displaystyle \frac{C_{NC2}^{-1}}{\frac{\overline{d}\frac{\sqrt{N\left(t\right)}}{\sqrt{\pi }}}{\overline{d}}-1}}=-{\displaystyle \frac{C_{NC2}^{-1}}{\frac{\sqrt{N\left(t\right)}}{\sqrt{\pi }}-1}}.$$

$C_{NC2}$ subsumes the effect of other factors in this rate, such as how the robots move and variations in the speed, similar to $C_{NC1}$ above. Hence,

$$\left({\displaystyle \frac{\sqrt{N\left(t\right)}}{\sqrt{\pi }}}-1\right)dN\left(t\right)\approx -C_{NC2}^{-1}dt\iff \int \left({\displaystyle \frac{\sqrt{N\left(t\right)}}{\sqrt{\pi }}}-1\right)dN\left(t\right)\approx -\int C_{NC2}^{-1}dt\iff$$

$${\displaystyle \frac{N{\left(t\right)}^{1.5}}{1.5\sqrt{\pi }}}-N\left(t\right)+n_0\approx -C_{NC2}^{-1}t,$$

for a constant $n_0$. As for $t=0$, $N\left(0\right)=N$, $n_0=-{\displaystyle \frac{N^{1.5}}{1.5\sqrt{\pi }}}+N,$ implying that when $N\left(t\right)=0$,

$$t\approx C_{NC2}\left({\displaystyle \frac{N^{1.5}}{1.5\sqrt{\pi }}}-N\right).$$

The desired result is obtained by summing parts (a)-(e) and simplifying the final result. □

3.2. Single Queue Former Algorithm

Algorithm Description: The Single Queue Former (SQF) algorithm [23] makes one queue shaped as a rectangular corridor that goes towards the target. This corridor has a width equal to the circular target diameter, $2s$, and a length equal to the working radius, D (Figure 6). The robots are permitted to enter the target region only by this queue. Potential fields are deployed to form this queue (Figure 6(a)) and to guide the robots to the target area exit (Figure 6(b)) efficiently. The corridor starts from the current target centre $\mathbf{G}=(G_x,G_y)$. Without any loss of generality, the corridor has vertices located at $(G_x+s,G_y+D)$, $(G_x-s,G_y+D)$, $(G_x-s,G_y)$ and $(G_x+s,G_y)$.

A rotational force field is applied to robots in the working area to enter the common target area through the corridor. The field rotation centre is at the target centre. The robots are submitted to an attractive force towards the target once they reach the corridor. After they arrive at the target area, another rotational force field is applied to them, whose centre is either at the leftmost point of the working area, $\mathbf{P}=(P_x,P_y)=(G_x-D,G_y)$, or at the rightmost point, $\mathbf{Q}=(Q_x,Q_y)=(G_x+D,G_y)$, depending on the position of the next target.

To control the attractive forces, each robot has the states going to the target ($GT$), leaving the target ($LT$) and going to the corridor ($GC$), which respectively means the robot is going straightly to the target region, it is leaving it, and it is going to the corridor. The starting state of the robots is $GT$ since they start outside the working area. Consider the conditions named W, O and A, represented as well-formed formulas:

$$\begin{array}{c}W=\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{G}\parallel \le D,\\ O=r_{iy}\left(t\right)<G_y\vee |r_{ix}\left(t\right)-G_x|>s,\\ A=\exists t^{\prime }\le t:\parallel {\mathbf{r}}_i\left(t^{\prime }\right)-\mathbf{G}\parallel \le s,\end{array}$$

where ${\mathbf{r}}_i\left(t\right)=(r_{ix}\left(t\right),r_{iy}\left(t\right))$ is the position of the robot i at time t. In other terms, W is true when the robot is inside the working area, O, when it is outside the corridor unlimited from above and A, when it has already arrived at the target area. Note that W and O do not depend on past time as A. Thus, the states of the robot can be represented as formulas using these conditions, which are mutually exclusive:

$$\begin{array}{c}GT=\neg W\vee \neg O\wedge \neg A,\\ GC=W\wedge O\wedge \neg A,\\ LT=W\wedge A.\end{array}$$

By using these conditions, the attractive force is represented as follows. Robots outside the corridor follow a force according to Passos et al. [23]

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)=\left\{\begin{array}{cc}{K}_a{\displaystyle \frac{(r_{iy}\left(t\right)-G_y,-r_{ix}\left(t\right)+G_x)}{\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{G}\parallel }},\hfill & \mathrm{if}GC\wedge r_{ix}\left(t\right)-G_x\le 0,\hfill \\ K_a{\displaystyle \frac{(-r_{iy}\left(t\right)+G_y,r_{ix}\left(t\right)-G_x)}{\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{G}\parallel }},\hfill & \mathrm{if}GC\wedge r_{ix}\left(t\right)-G_x>0,\hfill \end{array}\right.$$

(4)

where $K_a$ is a constant for setting the force magnitude. In other words, for the left-hand side of the circular working area, a clockwise rotational field is applied, and for the right-hand side, an anti-clockwise one.

However, if the robots are in the corridor, they go towards the target area following the same attractive force equation as in Equation 2, that is,

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)={\displaystyle \frac{K_a(\mathbf{G}-{\mathbf{r}}_i\left(t\right))}{\parallel \mathbf{G}-{\mathbf{r}}_i\left(t\right)\parallel }},\mathrm{if}GT.$$

(5)

Robots exiting nearby the target are constrained by another rotational field. This field also is only applied inside the circular working area:

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)=\left\{\begin{array}{cc}{K}_a{\displaystyle \frac{(r_{iy}\left(t\right)-P_y,-r_{ix}\left(t\right)+P_x)}{\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{P}\parallel }},\hfill & \mathrm{if}LT\wedge G_x^{\prime }\le G_x,\hfill \\ K_a{\displaystyle \frac{(-r_{iy}\left(t\right)+Q_y,r_{ix}\left(t\right)-Q_x)}{\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{Q}\parallel }},\hfill & \mathrm{if}LT\wedge G_x^{\prime }>G_x,\hfill \end{array}\right.$$

(6)

where ${\mathbf{G}}^{\prime }=(G_x^{\prime },G_y^{\prime })$ is the new target location. In other terms, for a robot with a new target located on the left-hand side of the previous target centre, a clockwise rotational field is applied and, on the right-hand side, an anti-clockwise one.

Finally, in addition to the default influence radius $I_d$, another constant is used by the robots when calculating repulsive forces: $I_{min}$, with $I_{min}<I_d$, the minimum influence radius allowed. For robots inside the corridor or exiting the target region, the influence radius is $I_{min}$. Now consider a robot outside the corridor but on the upper half of the circular working area. Let $d^{\prime }$ be the distance between the robot and a vertical line in the middle of the corridor (Figure 6(b)). Its influence radius I varies in relation to $d^{\prime }$ and is set to $I=I_{min}+d^{\prime }$ only for $0\le d^{\prime }\le I_d-I_{min}$. This range for $d^{\prime }$ guarantees that $I_{min}\le I\le I_d$. For the other robots, the influence radius is $I_d$. Thus, the influence radius and repulsive force are given by

$$\begin{array}{c}I\left({\mathbf{r}}_i\left(t\right)\right)=\left\{\begin{array}{cc}{I}_{min},\hfill & \mathrm{if}W\wedge (GT\vee LT),\hfill \\ I_{min}+|r_{ix}\left(t\right)-G_x|,\hfill & \mathrm{if}GC\wedge r_{iy}\left(t\right)>G_y\wedge |r_{ix}\left(t\right)-G_x|\hfill \\ \hfill & \phantom{\mathrm{if}}<I_d-I_{min},\hfill \\ I_d,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\\ {\mathbf{F}}_{\mathbf{R}}({\mathbf{r}}_i\left(t\right),{\mathbf{r}}_j\left(t\right))=\left\{\begin{array}{cc}{K}_r\left({\displaystyle \frac{1}{\parallel {\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j\left(t\right)\parallel }}-{\displaystyle \frac{1}{I\left({\mathbf{r}}_i\left(t\right)\right)}}\right){\displaystyle \frac{{\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j\left(t\right)}{\parallel {\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j{\left(t\right)\parallel }^3}},\hfill & \mathrm{if}\parallel {\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_j\left(t\right)\parallel \hfill \\ \hfill & \phantom{\mathrm{if}}\le I\left({\mathbf{r}}_i\left(t\right)\right),\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

(7)

Observe that plugging those attractive and repulsive forces in Equations 4–7 into the movement controller (Equation 1) yields an intricate differential equation with cases that not only depend on the position of an individual robot for attraction but also the robots in the neighbourhood due to the repulsive force and past time. Moreover, the influence radius calculated in Equation 7 also has different cases, expanding the chain of related expressions to calculate the solutions of the differential equations for each robot. Due to this complexity, the inference of the SQF algorithm task completion time function from the controller equations becomes more difficult than when the robots use NC. As a result, the estimated task completion time function equation is calculated by observing the behaviour of the robots in experiments, similarly to the previous section.

Samples of the Observations: As an example, Figure 7 and Figure 8 illustrate an experiment of the SQF algorithm with the default value of its parameter $I_{min}=1$ m. Red, green and yellow robots are in the state going to the target, going to the corridor and leaving the target, respectively. Black robots have exited the working area and are going to the next target. The starting positions are shown in Figure 7(a). The robots go in the direction of the target and change to state going to the corridor as they reach the working area (Figure 7(b)). The longest time for all the robots to enter the corridor is approximately the time for the bottom-most robot to go to it over the working area border (marked with a circle and an arrow in Figure 7(c). Its time includes the waiting of all robots ahead due to those occupying the corridor while they move towards the target area. After that, all robots can access the corridor (Figure 7(d)). After the last robot reaches the target area (marked with a circle and an arrow in Figure 8(e)), all robots follow the exit force field and are on different sides depending on their next target (Figure 8(f)).

Depending on the target size, the area for leaving the target may have robots still going to the entrance corridor. In this circumstance, the robots going to the target area must wait before they reach a location without any robot leaving the working area because these robots keep passing through them. Figure 9 and Figure 10 present an example of such a situation. Figure 9(a) illustrates the initial configuration. Figure 9(b) presents the first robot to reach the target area through the corridor (marked with a circle and an arrow). Figure 9(c) shows a robot that has arrived in the target area without going through the corridor and is leaving the target area. As the space between the target area and the working area border is small, the other robots push the robots near the target, and some may go to the target area without passing through the corridor. By the time the first robot leaves the working area (marked with a circle and an arrow in Figure 9(d)), there are robots going to the corridor in green and robots reaching the working area in red in the middle of the way induced by the potential field to leave the target area. These robots will take more time to go to the corridor than when s was smaller. Thus, they create a barrier reducing the area for the robots to leave, and the congestion also blocks the robots going to the corridor, thus, increasing their time to leave (Figure 10(e)). Figure 10(f) shows the end of this experiment. Observe that most of the robots are not in this figure. The robots on it are examples of robots that took more time to go to the corridor because they were on the leaving route when the other robots were leaving, as described above. Hence, these robots in Figure 10(f) were the last to overcome the robots exiting through the leaving area.

As before, an approximation of the expected task completion time for the SQF and its explanation are presented as follows.

Estimation 4.

The estimated expected time for N robots starting at a random position with a distance from the target centre in $[D,D+E$ to arrive at the common target area and leave the working area using the SQF algorithm is

$$\begin{array}{cc}\hfill \mathbb{E}\left[t_{SQF}\right]& \approx {\displaystyle \frac{\sqrt{D^2-s^2}}{\overline{v}}}\left(arcsin\left({\displaystyle \frac{s}{D}}\right)+arccos\left({\displaystyle \frac{\sqrt{D^2-s^2}}{2D}}\right)\right)+{\displaystyle \frac{D-s}{\overline{v}}}\hfill \\ \hfill & +{\displaystyle \frac{N(\pi -arcsin(s/D\left)\right)D}{(N+2)\overline{v}}}+{\displaystyle \frac{NE}{\overline{v}(N+1)}}+wai{t}_2\left(N\right),\hfill \end{array}$$

(8)

for

$$wai{t}_2\left(N\right)=\left\{\begin{array}{cc}{C}_{SQF1}N,\hfill & \mathrm{if}N\le {\displaystyle \frac{1}{2{\overline{d}}^2}}\left(2D^{2}arccos\left({\displaystyle \frac{s}{D}}\right)+2s^2+2{(D-s)}^{2}arccos\left({\displaystyle \frac{D-s}{2D}}\right)\right.\hfill \\ & \left.\phantom{if}+4D^{2}arcsin\left({\displaystyle \frac{D-s}{2D}}\right)-(D-s)\sqrt{4D^2-{(D-s)}^2}\right),\hfill \\ C_{SQF2}{N}^2,\hfill & \mathit{otherwise}\hfill \end{array}\right.$$

and constants $C_{SQF1}$ and $C_{SQF2}$.

Explanation. 

The estimation of the expected time to complete the SQF algorithm is better explained by analysing the experiment backwardly from the last robot to complete the task. Its calculation has five parts: the expected time for (a) the last robot leaving the working area from the target area (e.g., the time from Figure 8(e) to 8(f)), (b) the last robot going from the topmost edge of the corridor to the target area (e.g., from Figure 7(d) to 8(e)), (c) robots waiting to access the corridor (e.g., Figures 7(c), 9(d) and 10(e)), (d) the bottom-most robot travelling from where it first entered the working area border to the corridor (e.g., the time from Figure 7(b) to Figure 7(d) but excluding the waiting time (c)), (e) the bottom-most robot going from its starting position to the working area (e.g., Figure 7(a)–7(b)).

Part (a) is the time from the expected location in the target area reached by the last robot to the border of the working area by following a circular arc due to the exit force field. On average, the starting position of this arc is at the point $A_1$ in the circular target area and ends at $A_2$ in Figure 11(a), which shows this arc in dashed line. This arc belongs to a circle with radius $\sqrt{D^2-s^2}$ as illustrated by the right triangle formed by the points $A_1$ (at the target area circle with radius s by straightly following the corridor corner $A_4$ in Figure 11(b)), the target centre located at $\mathbf{G}$ and the point at $\mathbf{P}$1. Its length is proportional to the angle ${\alpha }_1+{\alpha }_2$, with ${\alpha }_1=arcsin(s/D)$ in the right triangle and ${\alpha }_2=arccos(\sqrt{D^2-s^2}/\left(2D\right))$ in the isosceles triangle $\mathbf{G}\mathbf{P}{A}_2$. Thus, the expected time (a) is given by $\sqrt{D^2-s^2}(arcsin(s/D)+arccos(\sqrt{D^2-s^2}/\left(2D\right)))/\overline{v}$.

Part (b) is given by the time to move from the topmost edge (see point $A_4$ in Figure 11(b)) to the target area, that is, $(D-s)/\overline{v}$.

Part (c) depends on the target area size and the number of robots. If the number of robots is greater than the entry area capacity, there will be robots going to the corridor in the leaving area while others are leaving the working area, causing congestion (as in Figure 9(d)–10(f)). The entry area capacity is calculated from the part of the working area likely to be occupied by robots going to the target area and unused for leaving. Thus, it includes the part of the upper half working area that the robots are probable to occupy inside the corridor, including half of the target area, but excluding the area used only by the robots that first arrived in the corridor (for instance, in Figure 9(d) and 10(e)). The entry area also includes two lateral areas where robots leaving the target area are unlikely to go, given the SQF force field for leaving from the left- and rightmost points of the target area. The entry area is shaded in Figure 11(c). For the calculation of the entry area, three areas are considered – namely, $Are{a}_1$, $Are{a}_2$ and $Are{a}_3$ – defined in the following paragraphs. Due to symmetry, $Are{a}_2$ and $Are{a}_3$ are doubled for the final outcome.

$Are{a}_1$ is the working area upper half minus the area of $△D_1{D}_2{D}_3$ and the area of the circular segment of points $D_1$ and $D_2$. For the area of $△D_1{D}_2{D}_3$, the height $h_1$ is given by the height of the isosceles triangle $\mathbf{G}{D}_1{D}_2$ – which is $\sqrt{D^2-s^2}$ from its inner right triangle – minus the radius s. Thus, $h_1=\sqrt{D^2-s^2}-s$, and the area of $△D_1{D}_2{D}_3$ is $s(\sqrt{D^2-s^2}-s)$. For the circular segment area, ${\beta }_1$ is obtained from the isosceles triangle $\mathbf{G}{D}_1{D}_2$ (with equal sides measuring D and base side, $2s$). Thus, ${\beta }_1=2arcsin(s/D)$, and the circular segment area is obtained by subtracting the area of $△\mathbf{G}{D}_1{D}_2$ from the working area sector of ${\beta }_1$: ${\displaystyle \frac{1}{2}}{D}^2{\beta }_1-{\displaystyle \frac{1}{2}}2s\sqrt{D^2-s^2}=D^{2}arcsin(s/D)-s\sqrt{D^2-s^2}$. Hence, $Are{a}_1={\displaystyle \frac{\pi }{2}}{D}^2-s(\sqrt{D^2-s^2}-s)-D^{2}arcsin(s/D)+s\sqrt{D^2-s^2}=D^2\left({\displaystyle \frac{\pi }{2}}-arcsin\left({\displaystyle \frac{s}{D}}\right)\right)+s^2=D^{2}arccos\left({\displaystyle \frac{s}{D}}\right)+s^2$.

$Are{a}_2$ and $Are{a}_3$ comprise the area on the left-hand side where robots leaving the target area are almost unlikely to pass, given the leaving force field from the leftmost point of the target area. $Are{a}_2$ is a circular sector area of angle ${\beta }_2$, and $Are{a}_3$ is the circular segment area of points $D_4$ and $D_5$ with angle ${\beta }_3$. ${\beta }_2$ and ${\beta }_3$ are obtained by the isosceles triangle $\mathbf{G}{D}_4{D}_5$ of equal sides measuring D and base side, $D-s$, obtained from the radius of the circular sector going from the leftmost point of the circular target area to the point $D_5$. Thus, ${\beta }_2=arccos\left({\displaystyle \frac{D-s}{2D}}\right)$ and ${\beta }_3=2arcsin\left({\displaystyle \frac{D-s}{2D}}\right)$. Consequently, $Are{a}_2={\displaystyle \frac{1}{2}}{(D-s)}^2{\beta }_2={\displaystyle \frac{1}{2}}{(D-s)}^{2}arccos\left({\displaystyle \frac{D-s}{2D}}\right)$, and $Are{a}_3$ is the circular sector area of angle ${\beta }_3$ minus the area of the isosceles triangle $\mathbf{G}{D}_4{D}_5$, that is, $Are{a}_3=D^{2}arcsin\left({\displaystyle \frac{D-s}{2D}}\right)-{\displaystyle \frac{D-s}{4}}\sqrt{4D^2-{(D-s)}^2}$. The entry capacity is approximated by considering each robot occupying a squared area of ${\overline{d}}^2$ ${\mathrm{m}}^2$ inside the entry area: ${\displaystyle \frac{1}{{\overline{d}}^2}}(Are{a}_1+2(Are{a}_2+Are{a}_3))={\displaystyle \frac{1}{2{\overline{d}}^2}}\left(2D^{2}arccos({\displaystyle \frac{s}{D}}\right)+2s^2+2{(D-s)}^{2}arccos\left({\displaystyle \frac{D-s}{2D}}\right)+4D^{2}arcsin\left({\displaystyle \frac{D-s}{2D}}\right)-(D-s)\sqrt{4D^2-{(D-s)}^2})$.

If the number of robots is less than or equal to the entry capacity, experiments show that the waiting time is linear in relation to the number of robots; otherwise, it is quadratic (due to the interference of the robots going to the corridor outside the entry area to the robots leaving the working area). Therefore, in the former case, part (c) is estimated by $C_{SQF1}N$, and, in the latter case, by $C_{SQF2}{N}^2$, for constants $C_{SQF1}$ and $C_{SQF2}$. As before, these constants abstract characteristics intrinsic to the robot dynamics and the environment where the algorithm is applied.

Part (d) is calculated from the arc length from the estimated initial position of the last robot to the corridor topmost edge (point $A_4$ in Figure 11(b)). In Figure 11(b), the point $A_3$ represents the maximum possible initial position. The estimated arc angle is the expected maximum from approximately $N/2$ robots on the left-hand side uniformly sampled over $0,\pi -{\alpha }_3]$, i.e., ${\displaystyle \frac{\frac{N}{2}(\pi -{\alpha }_3)}{\frac{N}{2}+1}}=$ ${\displaystyle \frac{N(\pi -{\alpha }_3)}{N+2}}$[68, p. 229]. From the right triangle formed by $A_4,A_5$ and $\mathbf{G}$ in Figure 11(b), ${\alpha }_3=arcsin(s/D)$. Thus, the estimated time of part (d) is $\left({\displaystyle \frac{N(\pi -arcsin(s/D\left)\right)}{N+2}}\right)D/\overline{v}$.

Part (e) is also obtained from the expected maximum but for N robots and the uniformly distributed distance from the starting point to the working area. It is given by ${\displaystyle \frac{NE}{\overline{v}(N+1)}}$. The desired result is obtained by summing parts (a)-(e). □

3.3. Touch and Run Vector Fields Algorithm

Algorithm Description: The Touch and Run Vector Fields (TRVF) algorithm [23] creates K lanes around the target area by using potential fields. Each lane goes from the working area to next to the target area through an entrance path in a straight fashion, makes a circular curve intersecting only one point of the target area through a circular curved path, and then leaves the working area directly through an exit path (Figure 12). An entrance path is parallel to the exit path on the opposing lane, and both are apart from d meters. Between them, a parallel ray passes from the target centre to the working area border. There are K such rays that divide the target area into sectors with angle $\alpha =2\pi /K$. Each lane is contained in such sectors called the central angle region.

Consider a central angle region $z=⌊{\displaystyle \frac{\eta }{\alpha }}⌋+1$, such that $\eta$ is the angle of the vector ${\mathbf{r}}_i\left(t\right)-\mathbf{G}$ with the x-axis. Each lane is formed by four waypoints: ${\mathbf{w}}_{\mathbf{1}}$ and ${\mathbf{w}}_{\mathbf{2}}$ forms the entrance path, ${\mathbf{w}}_{\mathbf{2}}$ and ${\mathbf{w}}_{\mathbf{3}}$, the circular curved path, and ${\mathbf{w}}_{\mathbf{3}}$ and ${\mathbf{w}}_{\mathbf{4}}$, the exit path (Figure 12).

Thus [23],

$${\mathbf{w}}_{\mathbf{1}}=\mathbf{G}+D(cos\left(z\alpha \right),sin\left(z\alpha \right))+{\displaystyle \frac{d}{2}}(sin\left(z\alpha \right),-cos\left(z\alpha \right));$$

$$\begin{array}{c}{\mathbf{w}}_{\mathbf{2}}=\mathbf{G}+\sqrt{{(r+s)}^2-{(r+d/2)}^2}(cos\left(z\alpha \right),sin\left(z\alpha \right))+\\ {\displaystyle \frac{d}{2}}(sin\left(z\alpha \right),-cos\left(z\alpha \right));\end{array}$$

$$\begin{array}{c}{\mathbf{w}}_{\mathbf{3}}=\mathbf{G}+\sqrt{{(r+s)}^2-{(r+d/2)}^2}(cos\left((z-1)\alpha \right),\\ sin\left((z-1)\alpha \right))+{\displaystyle \frac{d}{2}}(sin\left((z-1)\alpha \right),-cos\left((z-1)\alpha \right));\mathrm{and}\end{array}$$

$$\begin{array}{c}{\mathbf{w}}_{\mathbf{4}}=\mathbf{G}+D(cos\left((z-1)\alpha \right),sin\left((z-1)\alpha \right))+\\ {\displaystyle \frac{d}{2}}(sin\left((z-1)\alpha \right),-cos\left((z-1)\alpha \right))\end{array}$$

for the radius r of the circular path from ${\mathbf{w}}_{\mathbf{2}}$ to ${\mathbf{w}}_{\mathbf{3}}$. Due to Lemma 3 Passos et al. [22], $r={\displaystyle \frac{ssin(\alpha /2)-d/2}{1-sin(\alpha /2)}}$. This path has centre located at $\mathbf{c}=\mathbf{G}+(r+s)\left(cos\left(\left(z-{\displaystyle \frac{1}{2}}\right)\alpha \right)\right.,$$\left.sin\left(\left(z-{\displaystyle \frac{1}{2}}\right)\alpha \right)\right)$ [23].

The robot follows a straight line force field from the entrance and exit paths and an anti-clockwise orbit guided by the force field for the curved path. Let ${\mathbf{w}}_{\mathbf{i}}$ and ${\mathbf{w}}_{\mathbf{f}}$ be two arbitrary initial and final waypoints, and $\mathbf{C}$ and R any centre and radius to follow an orbit. The force field for the straight line following is expressed by (adapted from Nelson et al. [69])

$${\xi }_l=\left\{\begin{array}{cc}{\xi }_f-{\displaystyle \frac{\pi }{2}}\rho ,\hfill & \mathrm{if}ϵ>{\displaystyle \frac{d}{5}},\hfill \\ {\xi }_f-\left({\displaystyle \frac{k_s\frac{\pi }{2}v}{k_{\omega }{\left(\frac{d}{5}\right)}^{k_s}}}\right){ϵ}^{k_s-1}sin\left(\xi \right)-{\displaystyle \frac{\pi }{2}}{\left({\displaystyle \frac{ϵ}{\frac{d}{5}}}\right)}^{k_s},\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

$${\mathbf{F}}_{\mathbf{L}}({\mathbf{r}}_i\left(t\right),{\mathbf{w}}_{\mathbf{i}},{\mathbf{w}}_{\mathbf{f}})=\left\{\begin{array}{cc}{K}_a(cos\left({\xi }_l\right),sin\left({\xi }_l\right)),\hfill & \mathrm{if}{\displaystyle \frac{({\mathbf{r}}_i\left(t\right)-{\mathbf{w}}_{\mathbf{i}})·{\mathbf{w}}_{\mathbf{fi}}}{\parallel {\mathbf{w}}_{\mathbf{fi}}{\parallel }^2}}\ge 1,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(9)

for ${\mathbf{w}}_{\mathbf{fi}}={\mathbf{w}}_{\mathbf{f}}-{\mathbf{w}}_{\mathbf{i}}$, ${\xi }_f=atan2(w_{fi,y},w_{fi,x})$, $ϵ=∥{\mathbf{r}}_i\left(t\right)-{\mathbf{w}}_{\mathbf{i}}-{\displaystyle \frac{({\mathbf{r}}_i\left(t\right)-{\mathbf{w}}_{\mathbf{i}})·{\mathbf{w}}_{\mathbf{fi}}}{\parallel {\mathbf{w}}_{\mathbf{fi}}{\parallel }^2}}{\mathbf{w}}_{\mathbf{fi}}∥$, $\rho =sign({\mathbf{w}}_{\mathbf{fi}}\times ({\mathbf{r}}_i\left(t\right)-{\mathbf{w}}_{\mathbf{i}}))$, a constant for setting the force magnitude $K_a$, the current orientation of the robot $\xi$, the constant for proportional angular speed controller $k_{\omega }$, the constant for exponentiation in the vector field calculation $k_s>1$ and the maximum linear speed v.

The force field for the orbit following is given by (adapted from Nelson et al. [69])

$${\xi }_c=\left\{\begin{array}{cc}\gamma -{\displaystyle \frac{5\pi }{6}}+{\displaystyle \frac{v}{\parallel \mathbf{q}\parallel }}sin(\xi -\gamma ),\hfill & \mathrm{if}\parallel \mathbf{q}\parallel >2R,\hfill \\ \gamma -{\displaystyle \frac{\pi }{2}}-{\displaystyle \frac{\pi }{3}}{\left({\displaystyle \frac{\parallel \mathbf{q}\parallel -R}{R}}\right)}^{k_o}-{\displaystyle \frac{v}{k_{\omega }\parallel \mathbf{q}\parallel }}sin(\xi -\gamma )\hfill & \\ -{\displaystyle \frac{k_{o}v\pi }{3R^{k_o}{k}_{\omega }}}{(\parallel \mathbf{q}\parallel -R)}^{k_o-1}cos(\xi -\gamma ),\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

$${\mathbf{F}}_{\mathbf{O}}({\mathbf{r}}_i\left(t\right),\mathbf{C},R,{\mathbf{w}}_{\mathbf{f}})=\left\{\begin{array}{cc}{K}_a(cos({\displaystyle \frac{\pi }{2}}-{\xi }_c),sin({\displaystyle \frac{\pi }{2}}-{\xi }_c)),\hfill & \mathrm{if}\mathbf{q}\times ({\mathbf{w}}_{\mathbf{f}}-\mathbf{C})\le 0,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(10)

for $\mathbf{q}=(q_x,q_y)={\mathbf{r}}_i\left(t\right)-\mathbf{C}$, $\gamma =atan2(q_x,q_y)$ and the constant for exponentiation in this vector field calculation $k_o>1$.

In addition, the robots have six states to indicate in which part of the lane they are: going to the target ($GT$), going to the entrance path ($GEP$), on the entrance path ($OEN$), entering through the curve path ($ECP$), leaving through the curve path ($LCP$) and on the exit path ($OEX$). The initial state is going to the target ($GT$).

Consider the conditions named W, A, B, $EN$ and $TA$, represented as well-formed formulas:

$$\begin{array}{c}W=\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{G}\parallel \le D,\\ A=\exists t^{\prime }\le t:\parallel {\mathbf{r}}_i\left(t^{\prime }\right)-\mathbf{G}\parallel \le s,\\ B=\exists t^{\prime }\le t:({\mathbf{r}}_i\left(t^{\prime }\right)-\mathbf{G})\times ({\mathbf{w}}_{\mathbf{1}}-\mathbf{G})\le 0,\\ EN=\exists t^{\prime }\le t:{\displaystyle \frac{({\mathbf{r}}_i\left(t^{\prime }\right)-{\mathbf{w}}_{\mathbf{1}})·({\mathbf{w}}_{\mathbf{2}}-{\mathbf{w}}_{\mathbf{1}})}{\parallel {\mathbf{w}}_{\mathbf{2}}-{\mathbf{w}}_{\mathbf{1}}{\parallel }^2}}\ge 1,\\ TA=\exists t^{\prime }\ge t_A:({\mathbf{r}}_i\left(t^{\prime }\right)-\mathbf{c})\times ({\mathbf{w}}_{\mathbf{3}}-\mathbf{c})\le 0.\end{array}$$

W and A means the same as in SQF in the previous section. B, $EN$ and $TA$ respectively means in some moment the robot left the entering orbit towards the entrance path, left the entrance path and left the curved path orbit next to target area. $t_A$ is the time that robot i arrived at target area, i.e., $t_A={min}_t\parallel {\mathbf{r}}_i\left(t\right)-\mathbf{G}\parallel \le s$.

As in the previous section, the states can be expressed as mutually exclusive formulas using these conditions:

$$\begin{array}{c}GT=\neg W\\ GEP=W\wedge \neg A\wedge \neg B\wedge \neg EN\wedge \neg TA,\\ OEN=W\wedge \neg A\wedge B\wedge \neg EN\wedge \neg TA,\\ ECP=W\wedge \neg A\wedge B\wedge EN\wedge \neg TA,\\ LCP=W\wedge A\wedge B\wedge EN\wedge \neg TA,\\ OEX=W\wedge A\wedge B\wedge EN\wedge TA.\end{array}$$

Using them and Equations 9 and 10 the attractive forces of the TRVF algorithm are expressed as follows. Robots in $GT$ are attracted to the current target by a force ${\mathbf{F}}_{\mathbf{A}\mathbf{1}}$ but are repelled by a ${\mathbf{F}}_{\mathbf{TR}}$ from the working area of the previous target position ${\mathbf{G}}^{\prime }$ (zero, if there is no previous target). The sum of these forces is normalised so that the resultant force has magnitude $K_a$:

$${\mathbf{F}}_{\mathbf{A}\mathbf{1}}\left({\mathbf{r}}_i\left(t\right)\right)={\displaystyle \frac{K_a(\mathbf{G}-{\mathbf{r}}_i\left(t\right))}{\parallel \mathbf{G}-{\mathbf{r}}_i\left(t\right)\parallel }},$$

$${\mathbf{F}}_{\mathbf{TR}}\left({\mathbf{r}}_i\left(t\right)\right)=\left\{\begin{array}{cc}{K}_r\left({\displaystyle \frac{1}{\parallel {\mathbf{G}}^{\prime }-{\mathbf{r}}_i\left(t\right)\parallel }}-{\displaystyle \frac{1}{D}}\right){\displaystyle \frac{{\mathbf{G}}^{\prime }-{\mathbf{r}}_i\left(t\right)}{\parallel {\mathbf{G}}^{\prime }-{\mathbf{r}}_i{\left(t\right)\parallel }^3}},\hfill & \mathrm{if}\parallel {\mathbf{G}}^{\prime }-{\mathbf{r}}_i\left(t\right)\parallel -D\le D,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)=K_a{\displaystyle \frac{{\mathbf{F}}_{\mathbf{A}\mathbf{1}}\left({\mathbf{r}}_i\left(t\right)\right)+{\mathbf{F}}_{\mathbf{TR}}\left({\mathbf{r}}_i\left(t\right)\right)}{\parallel {\mathbf{F}}_{\mathbf{A}\mathbf{1}}\left({\mathbf{r}}_i\left(t\right)\right)+{\mathbf{F}}_{\mathbf{TR}}\left({\mathbf{r}}_i\left(t\right)\right)\parallel }},\mathrm{if}GT.$$

(11)

The robots in $OEN$ and $OEX$ only follow a straight line force field in ${\mathbf{F}}_{\mathbf{L}}$ with different parameters:

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)=\left\{\begin{array}{cc}{\mathbf{F}}_{\mathbf{L}}({\mathbf{r}}_i\left(t\right),{\mathbf{w}}_{\mathbf{1}},{\mathbf{w}}_{\mathbf{2}}),\hfill & \mathrm{if}OEN,\hfill \\ {\mathbf{F}}_{\mathbf{L}}({\mathbf{r}}_i\left(t\right),{\mathbf{w}}_{\mathbf{3}},{\mathbf{w}}_{\mathbf{4}}),\hfill & \mathrm{if}OEX.\hfill \end{array}\right.$$

(12)

Robots in $GEP$ are only guided by an orbital force field ${\mathbf{F}}_{\mathbf{O}}$ with the following parameters:

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)={\mathbf{F}}_{\mathbf{O}}({\mathbf{r}}_i\left(t\right),\mathbf{G},D,{\mathbf{w}}_{\mathbf{1}}),\mathrm{if}GEP.$$

(13)

In $ECP$ and $LCP$, they are submitted to a sum of an orbital force and a stronger attractive force, and the resultant modulus of this sum is constrained to $K_a$:

$${\mathbf{F}}_{\mathbf{A}\mathbf{2}}\left({\mathbf{r}}_i\left(t\right)\right)=1.5{\displaystyle \frac{K_a(\mathbf{G}-{\mathbf{r}}_i\left(t\right))}{\parallel \mathbf{G}-{\mathbf{r}}_i\left(t\right)\parallel }},$$

$${\mathbf{F}}_{\mathbf{A}\mathbf{3}}\left({\mathbf{r}}_i\left(t\right)\right)=1.5{\displaystyle \frac{K_a({\mathbf{w}}_{\mathbf{3}}-{\mathbf{r}}_i\left(t\right))}{\parallel {\mathbf{w}}_{\mathbf{3}}-{\mathbf{r}}_i\left(t\right)\parallel }},$$

$${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)=\left\{\begin{array}{cc}{K}_a{\displaystyle \frac{{\mathbf{F}}_{\mathbf{A}\mathbf{2}}\left({\mathbf{r}}_i\left(t\right)\right)+{\mathbf{F}}_{\mathbf{O}}({\mathbf{r}}_i\left(t\right),\mathbf{c},r,{\mathbf{w}}_{\mathbf{3}})}{\parallel {\mathbf{F}}_{\mathbf{A}\mathbf{2}}\left({\mathbf{r}}_i\left(t\right)\right)+{\mathbf{F}}_{\mathbf{O}}({\mathbf{r}}_i\left(t\right),\mathbf{c},r,{\mathbf{w}}_{\mathbf{3}})\parallel }},\hfill & \mathrm{if}ECP,\hfill \\ K_a{\displaystyle \frac{{\mathbf{F}}_{\mathbf{A}\mathbf{3}}\left({\mathbf{r}}_i\left(t\right)\right)+{\mathbf{F}}_{\mathbf{O}}({\mathbf{r}}_i\left(t\right),\mathbf{c},r,{\mathbf{w}}_{\mathbf{3}})}{\parallel {\mathbf{F}}_{\mathbf{A}\mathbf{3}}\left({\mathbf{r}}_i\left(t\right)\right)+{\mathbf{F}}_{\mathbf{O}}({\mathbf{r}}_i\left(t\right),\mathbf{c},r,{\mathbf{w}}_{\mathbf{3}})\parallel }},\hfill & \mathrm{if}LCP.\hfill \end{array}\right.$$

(14)

Again, the controller in Equation 1 with the attractive force described by Equations 11–14 results in a complicated differential equation with mutual dependence. Consequently, the estimated expected task completion time function equation for the TRVF algorithm is also calculated by observing the experiments.

Sample of the Observations: As an example, Figure 13 and Figure 14 present an execution of the TRVF algorithm with its default values (Table 3). Robots in red, cyan, blue, magenta, yellow, orange and black represent the states $GT$, $GEP$, $OEN$, $ECP$, $LCP$, $OEX$ and that they left the previous target area, respectively. Figure 13(a) illustrates the starting configuration. Robots go to the entrance path conducted by an orbital force around the working area. Then, they proceed to the target through the entrance path (Figure 13(b)). Figure 13(c) shows the first robot to reach the target in this experiment (marked with a circle and an arrow). Depending on the number of robots, some robots still have to wait to enter the entrance path (in cyan) while the first one to leave the working area exited (the black robot with a circle and an arrow in the bottom of Figure 13(d)). After that wait, the last robot to access the working area proceeds through the entrance path in the lane, arrives in the target area (with a circle and an arrow in Figure 14(e)) and leaves by the exit path, facing a few robots going in different directions (Figure 14(f)).

From such examples, an approximation of the expected task completion time for the TRVF and its explanation are shown below.

Estimation 5.

The estimated expected time for N robots starting at a random position with a distance from the target centre in $[D,D+E$ to arrive at the common target area and leave the working area using the TRVF algorithm is

$$\mathbb{E}\left[t_{TRVF}\right]\approx {\displaystyle \frac{E}{\left(\frac{arctan\left(\frac{d}{2D}\right)}{2\pi }N+1\right)\overline{v}}}+{\displaystyle \frac{2d_s+r(\pi -\alpha )}{\overline{v}}}+C_{TRVF}{\displaystyle \frac{N}{K}},$$

(15)

for a constant $C_{TRVF}$, $\alpha =2\pi /K$, $r={\displaystyle \frac{ssin(\alpha /2)-d/2}{1-sin(\alpha /2)}}$ and $d_s=\sqrt{D^2-d^2/4}-\sqrt{s(2r+s)-d\left(r+{\displaystyle \frac{d}{4}}\right)}$.

Explanation. 

To estimate the expected time to complete the task with the TRVF algorithm, three parts are needed: the expected time for (a) the first robot reaching an entrance path without moving across the border of the working area circular sector from its starting point, (b) the first robot going to the target area and then leaving it until reaching the working area border by one of the K lanes (e.g., Figure 13(b)–Figure 13(d)), (c) the waiting and travel of the last robot through the lane.

Part (a) is approximated by the expected minimum distance of a robot from its starting position going directly to an entrance path of a lane without moving next to the border of the working area circular sector. This approximation considers that it is more probable that the first robot to arrive at the entrance path was initially in front of the entrance path. Otherwise, the robot would have to move next to the working area border until it arrives at the entrance path, which is more distant, as the cyan robots in Figure 13(b)–13(d). From Figure 15(a), the sector angle in front of the entrance path is given by ${\alpha }_4=arctan(d/\left(2D\right))$. Thus, the expected time is calculated from the expected minimum distance [68, p. 229] for the number of robots located in the sector angle ${\alpha }_4$ when the experiment started (${\displaystyle \frac{{\alpha }_4}{\alpha }}{\displaystyle \frac{N}{K}}$), i.e.,

$${\displaystyle \frac{E}{\left(\frac{{\alpha }_4}{\alpha }\frac{N}{K}+1\right)\overline{v}}}={\displaystyle \frac{E}{\left(\frac{{\alpha }_4}{2\pi }N+1\right)\overline{v}}}={\displaystyle \frac{E}{\left(\frac{arctan\left(\frac{d}{2D}\right)}{2\pi }N+1\right)\overline{v}}}.$$

Figure 15(b) helps to understand how part (b) is calculated. This part is obtained from the time of the first robot moving at average speed through a distance of $d_s$ meters via the entrance path, going through the circular path until the point $A_6$, leaving by this curved path and running again a distance of $d_s$ meters via the exit path. The curved path length is obtained from the circular sector angle $\beta =\pi -\alpha$ (because the sum of the angles inside the quadrilateral containing the point $A_6$ must be $2\pi$ and the internal angle on the right-hand side is equal to $\alpha$ due to parallelism of the line in the lane).

The length $d_s$ in Figure 15(b) is calculated from $\overline{C_{1}C}$ and $e_s=\left|\overline{A_1{C}_1}\right|$ in Figure 15(c). The triangle $ABC$ is the same depicted in the proof in [22, Lemma 2], where the distance to the target centre for the robot to start turning was calculated using this triangle, so, from the proof of that lemma,

$$|\overline{AC}|=\sqrt{{(r+s)}^2-{(r+d/2)}^2}=\sqrt{s(2r+s)-d\left(r+{\displaystyle \frac{d}{4}}\right)}.$$

From the right triangle $A{B}_1{C}_1$, $|\overline{A{C}_1}|=\sqrt{D^2-d^2/4}$. As $D=|\overline{A_{1}A}|=|\overline{A{C}_1}|+|\overline{A_1{C}_1}|=\sqrt{D^2-d^2/4}+e_s$, $e_s=D-\sqrt{D^2-d^2/4}$. Also, $D=|\overline{A_{1}A}|=|\overline{A_1{C}_1}|+|\overline{C_{1}C}|+|\overline{AC}|=e_s+d_s+\sqrt{s(2r+s)-d\left(r+{\displaystyle \frac{d}{4}}\right)}$. Therefore,

$$d_s=D-e_s-\sqrt{s(2r+s)-d\left(r+{\displaystyle \frac{d}{4}}\right)}=\sqrt{D^2-d^2/4}-\sqrt{s(2r+s)-d\left(r+{\displaystyle \frac{d}{4}}\right)}.$$

As mentioned before, part (b) is the time to go through the entrance lane, the curved path and the exit path, measuring $d_s$, $r\beta$ and $d_s$, respectively. Consequently, part (b) is given by ${\displaystyle \frac{2d_s+r(\pi -\alpha )}{\overline{v}}}$.

From experiments, part (c) is proportional to the number of robots in each lane, that is, $C_{TRVF}{\displaystyle \frac{N}{K}}$, for a constant $C_{TRVF}$ that abstracts the characteristics of the experiment. The final result holds by adding parts (a)-(c). □

3.4. Mixed Teams

Algorithm Description: Mixed Teams (MT) occurs when multiple groups of swarms act in the same environment, but they do not know the algorithm used by other groups. Assume two groups of robots, group A of $M<N$ robots and group B of $N-M$ robots. As explained below, the estimation for MT shows that the estimations for the algorithms can be used as a basis for further derivations. A robot of group A at position ${\mathbf{r}}_i\left(t\right)$ executes an algorithm $Alg$ which returns a vector ${\mathbf{F}}_{Alg}$. Thus, their attractive force is simply ${\mathbf{F}}_{\mathbf{A}}\left({\mathbf{r}}_i\left(t\right)\right)={\mathbf{F}}_{Alg}\left({\mathbf{r}}_i\left(t\right)\right).$ As noticed in the previous algorithms, deriving an exact equation of the task completion time function from the controller equations is complicated because two groups of robots are executing possibly two different algorithms.

As shown in Appendix A, using NC has the best results among the experimented algorithms for the group A when the percentage of robots in group A to the total number of robots (${\displaystyle \frac{N}{M}}·100\%$) is lower than 30% and 60% for non-holonomic and holonomic robots, respectively. Thus, in the analysis below ${\mathbf{F}}_{Alg}$ is the same as in (2), and $Alg$ = NC was used in the following illustrations.

Samples of the Observations: Figure 16, Figure 17, Figure 18 and Figure 19 show the execution of MT by 10% and 50% of robots in group A (in grey) to the total number of robots, respectively, when the robots in group B are executing the SQF algorithm with its default values. In these figures, the total number of robots is $N=100$. In Figure 16(b), the usual behaviour of the 10% of the robots in group A is to go directly to the target. In Figure 16(c), some grey robots left the target region and are going to the next target on both sides (marked with circles and arrows). When the two last robots running NC leave the target area in Figure 16(d) (with circles and arrows), they try to go to the next target on the left-hand side, but other robots using the SQF are blocking their way. Due to that blockage, the two robots will continue through the SQF-induced leaving route until they can go to the left (Figure 17(e)). After that, only robots using SQF are in the experiment (Figure 17(f)). When few robots use NC, this often happens.

Figure 18 and Figure 19 have more robots using NC. From the bottom of Figure 18(b), more of them go directly to the target than in the previous example. Robots trying to reach the target may be pushed by robots using SQF on the way to leave the working area (Figure 18(c)). The cluttering formed by them may continue until the last robot using SQF reaches the target area (marked with a circle and an arrow in Figure 18(d)). This tendency to clutter occurs proportionally to the number of robots using NC as the algorithm of group A. The robots using NC which are oriented to the target but were on the SQF-induced leaving route are more likely to be the last ones to arrive at the target region (as the two grey robots marked with circles and arrows near the target in Figure 19(e)). At the end of this experiment, they are the last robots to leave the working area (one of them is marked with a circle and an arrow in Figure 19(f)).

Figure 20, Figure 21, Figure 22 and Figure 23 illustrate the execution of MT with 10% and 50% of robots in group A to the total number of robots for an experiment with the TRVF with the default parameter values as the algorithm of the robots in group B. In these figures, $N=100$. As occurred for SQF, using 10% of the robots has almost the same result as using only TRVF for all robots. Figure 20(b) shows the robots using NC arriving at the target area through the free space between the lanes. If a robot using NC is on the left-hand side, but its new target is on the right-hand side, it has to wait for robots using TRVF open space in the lane blocking its passage (as the grey robot marked with a circle and an arrow on the left side near the target area in Figure 20(c)). When a few robots use NC, they leave the working area faster than robots using TRVF, as they tend to follow the shortest distance to the target through the spaces between lanes (Figure 20(d)). In Figure 21(e), all robots executing NC have finished the task. Only robots with TRVF are leaving the target area by exit paths until the end of the experiment in Figure 21(f).

Figure 22 and Figure 23 exemplify the usage of 50% of robots in group A, and group B executes TRVF with default values. As the proportion of robots using NC is more than in the previous example, they are more likely to cover the open space near the target (Figure 22(b)) and cause congestion inside the target area (Figure 22(c)). If a lane has more robots using NC, they hamper robots using TRVF to go through it, as the exit lane on the top left-hand side in Figure 22(d). As for the 10% case, the last robot using NC leaves before the TRVF robots (marked with a circle and an arrow in Figure 23(e)). However, fewer TRVF robots than in 10% case are inside the working area until the end of the experiment (Figure 23(f)).

These illustrations hint that the proportion of robots in group A induces a cluttering similar to that occurred by using only NC. Thus, the task completion time tends to be near NC time as the number of robots in group A M increases and the algorithm followed by the robots in group B as M decreases. To better show this, Figure 24 and Figure 25 show the result from the experiments with MT for scenarios with the robots executing SQF and TRVF with different values of M. The bars mean the 99% confidence interval of the average for 40 runs for each value in the horizontal axis. (The examples above are shown for convenience for the explanation below. More detailed results are shown in Section 4.) Based on these observations, the next estimation shows an approximation for an MT, followed by its explanation.

Estimation 6.

The estimated expected time for N robots starting at a random position with a distance from the target centre in $[D,D+E$ to arrive at the common target area and leave the working area when $M<N$ robots are using MT with NC as the group A algorithm depends on the ratio $p=M/N$ and is expressed by

$$\begin{array}{cc}\hfill \mathbb{E}\left[t_{MT}\right]\left(p\right)& \approx C_{MTNC}\mathbb{E}\left[t_{NC}\right]+C_{MTB}\mathbb{E}\left[t_{AlgB}\right]\hfill \end{array}$$

(16)

for constants $C_{MTNC}$ and $C_{MTB}$, $AlgB$ being the algorithm followed by the $N-M$ robots in group B and its estimated expected time $\mathbb{E}\left[t_{AlgB}\right]$.

Proof 

(Explanation). Let $C_{MTNC}$ be a given constant abstracting the environment and the dynamics of the robots in group A for the ratio p and $C_{MTB}$ a similar constant for the robots in group B. Thus, as observed in the experiments (Figure 24 and Figure 25, Figure 35 and Figure 36 in Section 4 and Figure A9–Figure A14 in Appendix B), the estimated expected time is the sum of the estimated expected time of NC and the control algorithm followed by the robots in group B each multiplied by the constants $C_{MTNC}$ and $C_{MTB}$, respectively. □

4. Experiments and Results

Experiments were performed in the Stage simulator [36] with the algorithms in the previous section using the same task scenario and the default values used in the examples of the last section. Only the number of robots N is varied from 20 to 300 with steps of 20 and the robotic movement type (holonomic and non-holonomic, to show the adaptability of the estimations). For each N, 40 executions were run. In the figures of this section, the bar intervals are the 99% confidence interval of the mean value of the variable considered. “Experiments” and “Estimation” stand for these intervals of the experimental data mean and the equation estimated from these experiments, respectively.

Moreover, the constants for the estimated time equations were obtained with the least squares method from the mean value of the experimental data. It is important to heed that the presented estimations produce curves that can be fitted to the simulation results by tuning the coefficients of the expressions. Thus, using wrong equations with constants obtained from the least square method does not work appropriately. For example, the Estimation 3 has two cases depending on N. If one of the cases is ignored and the other is applied for every N, even by fitting its constant, the estimation will be further from the experiment data than the presented equation.2

Figure 26 and Figure 27 present the comparison between the experimental data and the estimated expected time computed by Estimation 3 for NC and $s\in \{3,5,7\}$ m. Observe that the estimation is closer as the value of s is smaller.

Figure 28 and Figure 29 show the comparison for the SQF algorithm with the estimated time in Estimation 4 for $s\in \{3,5,7\}$ m. The negative value of $C_{SQF1}$ in the non-holonomic case for $s=3$ m (Figure 28(b)) means that the waiting time was less than the predicted travel time of the last robot (i.e., the value from Estimation 4 without the term $C_{SQF1}N$). Notice that the shapes of the graphs are different for different type of robotic movement, although the equations are the same. This is because the average linear speed varies further with more robots in the non-holonomic case than with the holonomic one, as it will be shown later.

Figure 30 and Figure 31 show the comparison for TRVF algorithm with the estimated time in Estimation 5 for $s\in \{3,5,7\}$ m. Observe that using this estimation for holonomic robots and $s=7$ m is not as close as for non-holonomic ones (Figure 31). In the algorithms where the lower number of robots differs more from the estimation, the robots had little mutual influence as they were scattered. These estimations are not close to the experimental data for these values because they were deduced from the global behaviour affected by the influence of their surroundings assuming that the algorithms were applied with a large number of individuals.

Figure 32 gives the results of the root mean squared error normalised by the y-axis range (NRMSE) for the estimations in Figure 26, Figure 27, Figure 28, Figure 29, Figure 30 and Figure 31. From it, these estimations are considered close to the experimental data because the NRMSE values were lower than 0.03. In Figure 32, the lowest error for holonomic robots was for SQF with $s=5$ m while for non-holonomic case, it was NC with $s=3$ m. NC has the highest error for $s=7$ in both types of robotic movement because the graph seems linear instead of a curve.

In addition, this work assumes a fixed working area. That being so, notice that the influence of the number of robots has a limit because a higher number will clutter the environment no matter the algorithm. Thus, the parameters related to the environment area affect these functions in this sense.

Although it is not explicit from Estimations 3, 4 and 5 for the estimated expected time, $\overline{d}$ and $\overline{v}$ depend on N. Figure 33 presents the relation between the average speed and the number of robots for the algorithms for $s=3$ m obtained from the experiments. Note that the range for the holonomic case is smaller (around 0.25 m/s) than for the non-holonomic robots (from 0.05 to 0.25 m/s). The holonomic robot can maintain the linear speed in all directions, while the non-holonomic one needs to go forward or backward, depending on the situation, yielding positive and negative values for the linear speed. When more robots are next to each other, this forward and backward movement intensifies in the non-holonomic case. In the holonomic case, the slight rise in the speed as the number of robots increases is due to the velocity of a holonomic robot being calculated for every direction, as they can move freely in any direction. Thus, as more robots are added, they are deviating more and consequently, more speed is included in the average.

Figure 34 shows the relation between the mean distance between the robots and the number of robots obtained from experiments. For both types of robots, the trends in the mean distance are similar: as the number of robots grows, they have to move nigher to each other in the algorithms.

Figure 35 presents an instance of the estimations for MT with the proportion $p=10\%$ of robots in group A and SQF as the control algorithm of the robots in group B. For reference, the task completion time from the experiments of SQF and NC is shown. Appendix B contains the figures for p from 20% to 90% in increments of 10%. Figure 37(a) shows the results of the NRMSE for the estimations for these values of p. From the NRMSE values, they match better for non-holonomic robots from 70 to 90%.

Figure 35. Comparison of the estimated expected time and the experimental data for MT with NC as group A’s algorithm, group B is executing the SQF algorithm and $p=10\%$ for (a) holonomic robots ($C_{MTNC}=-0.0042$ and $C_{MTB}=0.9779$) and (b) non-holonomic robots ($C_{MTNC}=0.0069$ and $C_{MTB}=1.027$).

Figure 35. Comparison of the estimated expected time and the experimental data for MT with NC as group A’s algorithm, group B is executing the SQF algorithm and $p=10\%$ for (a) holonomic robots ($C_{MTNC}=-0.0042$ and $C_{MTB}=0.9779$) and (b) non-holonomic robots ($C_{MTNC}=0.0069$ and $C_{MTB}=1.027$).

Figure 36. Comparison of the estimated expected time and the experimental data for MT with NC as group A’s algorithm, group B is executing the TRVF algorithm and $p=10\%$ for (a) holonomic robots ($C_{MTNC}=0.0309$ and $C_{MTB}=0.9459$) and (b) non-holonomic robots ($C_{MTNC}=0.0766$ and $C_{MTB}=0.9323$).

Figure 36. Comparison of the estimated expected time and the experimental data for MT with NC as group A’s algorithm, group B is executing the TRVF algorithm and $p=10\%$ for (a) holonomic robots ($C_{MTNC}=0.0309$ and $C_{MTB}=0.9459$) and (b) non-holonomic robots ($C_{MTNC}=0.0766$ and $C_{MTB}=0.9323$).

Figure 37. Normalised root mean square error of the estimations with different percentages of robots in group A using NC for an MT with group B using (a) SQF and (b) TRVF.

Figure 37. Normalised root mean square error of the estimations with different percentages of robots in group A using NC for an MT with group B using (a) SQF and (b) TRVF.

Similarly, Figure 36 shows an example of the estimations for MT with TRVF as control algorithm of group B with $p=10\%$. As before, they have the task completion time from the experiments of TRVF and NC for reference. The figures for p ranging from 20% to 90% in steps of 10% appear in Appendix B. Figure 37(b) shows the results of the NRMSE for these estimations. Due to the estimation equation for TRVF working better for a higher number of robots, the ad hoc estimation reflects this result. Also, the NRMSE is higher for holonomic case.

Figure 38 displays the change in constants $C_{MTB}$ and $C_{MTNC}$ by the experimented percentage of robots in group A shown above for each tested control algorithm of the robots in group B and robots’ movement type. Note that when $C_{MTNC}$ rises, $C_{MTB}$ decreases because it is expected that these constants are proportional to the number of robots executing the algorithm associated with each of them. Although the values of the constants are different, the trend for the SQF algorithm is similar for both robots’ movement types (Figure 38(a)). The same happens for the TRVF algorithm in Figure 38(b).

5. Discussion

In this section, limitations, challenges and future improvements to the presented work are discussed. Firstly, this paper assumes that the only random operation is the initial position of the robots. Aside from that, the robots would behave deterministically. Therefore, the average task completion is due to the random initialisation. Thus, a best/worst-case scenario may be estimated by taking the standard deviation from above and below the average time in most cases. Although this is true for the simulations discussed in this paper, for real robots, there would be more stochastic events that should be inputted in these best-case/worst-case scenarios.

Secondly, the presented estimations are valid for changes in the parameters N and s as tested and described in Section 4. However, the estimations were not experimented with different values of every global parameter other than N and s (Table 2). The authors believe the estimations work for reasonable changes in the other global parameters. Nonetheless, significant changes were not tested and could cause unexpected results, since the underlying algorithms could have undesired behaviour with a wrong parametrisation.

Thirdly, the presented algorithms may be adapted to three-dimensional robots and experimented with their type of movements due to the recently increasing usage of aerial robots. Because this paper considered robots and the target regions in a plane, the equations must be extended to three dimensions. Three-dimensional robots have different types of movement that may influence the presented algorithms differently. For instance, aircraft-like robots cannot stop in mid-air, and some aerial robots are controlled by Dubins paths.

In addition, adapting a method made for a space with a fixed number of dimensions may not be effortless. The algorithms SQF and TRVF were born from a tentative by the authors to apply the queuing theory. However, this was unsuccessful because, in two-dimensional space, the queues behave differently depending on the path followed by the robots. In other words, the original queuing theory is one-dimensional, and adapting to two dimensions did not seem to be straightforward. Thus, only extending the algorithms may not be sufficient, but they also must be tested for feasibility. Should the extended algorithms be ineffective in three dimensions, new algorithms ought to be devised.

Another future step is to apply methods from statistical mechanics, in special for many-body problems, to this problem and to compare with those developed here. Such solutions may help to find a potential field that minimises the relation of the task completion time versus the number of robots. Additionally, these studies may answer whether finding such a potential field is possible. Notice that if this is possible, a time relation for a given field has to be given to be minimised, and this work focused on approximating such relation by exploring the usage of the geometry of the overall behaviour of the robots. These ideas are the basis for future improvements, which include finding a better potential field minimising the task completion time per number of robots, presuming that finding it is decidable.

Finally, the method of analysing the geometric behaviour from the experiments to approximate local-to-global functions used here may be applied to different problems of robotic swarms if a fixed number of robots follow similar paths on average for different starting positions. For instance, a patrolling swarm may go to specific areas regularly using a predefined route on average. By aggregating the routes, analysing their geometry and averaging by the frequency of usage, the expected time between patrols may be estimated. A similar analysis may be performed on a swarm covering an area to estimate the average time for completing that task. Another example is to approximate the volume per time of a swarm collecting waste by calculating the average time of collecting and accessing the areas with waste and the volume of waste from them by geometrically examining the used paths from the experiments and averaging by their frequency.

6. Conclusions

This paper presented estimations of the expected task completion time in relation to the number of robots in algorithms for the common target problem using potential fields. These relations were calculated by analysing the experiments with NC, SQF, TRVF and MT. Experiments with MT of two groups showed that if the percentage of robots in one group is lower than 30% and 60% for non-holonomic and holonomic robots, respectively, the results on the average total simulation time are better for NC as algorithm of a group, although it is not suitable for a swarm with only one group.

Besides the number of robots, the expected task completion time estimations are also expressed in terms of the average speed of the robots and algorithm parameters such as the working area radius, the target area radius, the number of lanes in TRVF, its central angle, the distance travelled through its lanes and the ratio of the number of robots in group A to the total number of individuals. In addition, the experimental relation of the average speed and distance between the robots as a function of the number of individuals in the swarm for SQF and TRVF was illustrated.

Therefore, this work is a fundamental first step to examining the common target problem as a local-to-global analysis of the task completion time from the specification of the individual robot to the swarm behaviour. An extensive and challenging study is needed to explain the local-to-global nature of these subsumed details. This invokes a theory to untangle this relationship, and this paper aims to start a discussion about that by presenting a baseline approximation to compare with future work. There is room for improvement on these estimations. Approaches using adaptations of the methods discussed in Section 2 (the perturbation theory, the coupled cluster technique, the density-functional theory, the mean-field theory and the variational method) are welcome to be compared.

Appendix A.1. Comparison of NC and SQF for the Group A and TRVF for Group B

Figure A1–Figure A3 illustrate the comparison of the total simulation time per total number of robots of an MT using SQF and NC as parameter $Alg$ with percentages of robots in group A to the total number of individuals ranging from 10% to 90% in steps of 10%. In these experiments, TRVF is the swarm control algorithm of group B. These figures are summarised in Figure A4, which shows the sum of the average total simulation time over each experimented number of robots by the percentage of robots in group A for MT. The bars correspond to the respective sum of the confidence interval ranges with $\rho =0.01$.

Figure A1. Comparison of the total simulation time by the number of robots for the experiments with MT using TRVF as the group B’s algorithm and SQF and NC as the group A’s algorithm by 10% of the robots for (a) holonomic and (b) non-holonomic robots, 20% for (c) holonomic and (d) non-holonomic robots and 30% for (e) holonomic and (f) non-holonomic robots.

Figure A1. Comparison of the total simulation time by the number of robots for the experiments with MT using TRVF as the group B’s algorithm and SQF and NC as the group A’s algorithm by 10% of the robots for (a) holonomic and (b) non-holonomic robots, 20% for (c) holonomic and (d) non-holonomic robots and 30% for (e) holonomic and (f) non-holonomic robots.

Figure A2. Comparison of the total simulation time by the number of robots for the experiments with MT using TRVF as the group B’s algorithm and SQF and NC as the group A’s algorithm by 40% of the robots for (a) holonomic and (b) non-holonomic robots, 50% for (c) holonomic and (d) non-holonomic robots and 60% for (e) holonomic and (f) non-holonomic robots.

Figure A2. Comparison of the total simulation time by the number of robots for the experiments with MT using TRVF as the group B’s algorithm and SQF and NC as the group A’s algorithm by 40% of the robots for (a) holonomic and (b) non-holonomic robots, 50% for (c) holonomic and (d) non-holonomic robots and 60% for (e) holonomic and (f) non-holonomic robots.

Figure A3. Comparison of the total simulation time by the number of robots for the experiments with MT using TRVF as the group B’s algorithm and SQF and NC as the group A’s algorithm by 70% of the robots for (a) holonomic and (b) non-holonomic robots, 80% for (c) holonomic and (d) non-holonomic robots and 90% for (e) holonomic and (f) non-holonomic robots.

Figure A3. Comparison of the total simulation time by the number of robots for the experiments with MT using TRVF as the group B’s algorithm and SQF and NC as the group A’s algorithm by 70% of the robots for (a) holonomic and (b) non-holonomic robots, 80% for (c) holonomic and (d) non-holonomic robots and 90% for (e) holonomic and (f) non-holonomic robots.

Figure A4. Comparison of the total simulation time sum over the number of robots by the percentage of robots in group A using NC and SQF as group A’s algorithms. The swarm control algorithm of the robots in group B is TRVF. (a) Holonomic robots. (b) Non-holonomic robots.

Figure A4. Comparison of the total simulation time sum over the number of robots by the percentage of robots in group A using NC and SQF as group A’s algorithms. The swarm control algorithm of the robots in group B is TRVF. (a) Holonomic robots. (b) Non-holonomic robots.

In Figure A1–Figure A3, for 30% to 70% for non-holonomic robots, NC as the algorithm of group A is significantly better on average until different values from 120 to 180 robots. Note that the number of robots where it starts to take less time decreases with the ratio of robots in group A.

Furthermore, SQF as the group A algorithm takes more or equal time on average than using NC for percentages from 10% to 20% of robots in group A using MT for both type of robots and 30% to 70% for holonomic robots with statistical significance (Figure A4). Although the SQF is faster, when used as the group A algorithm for a low percentage of robots, they follow a trajectory longer than they would if NC is in use. As fewer robots are going directly to the target region and there is less congestion, NC is better for low percentages. The variance is higher for SQF because the SQF robots leaving the target area may push the robots using TRVF away from the working area in the bottom lane for a long time.

Additionally, the usage of holonomic robots takes more time and varies more than non-holonomic robots for the SQF algorithm. Holonomic robots are faster to be repelled by another robot due to their lack of restrictions on moving in any direction. Because they can be repelled for longer periods, they tend to take more time going away from the working area.

Appendix A.2. Comparison of NC and TRVF for the Group A and SQF for Group B

Figure A5–Figure A7 show the comparison of NC and TRVF as the group A algorithms for ratios varying from 10% to 90% in steps of 10% for holonomic and non-holonomic robots. In these experiments, SQF is the swarm control algorithm of group B. As before, these figures are summarised in Figure A8, which presents the sum over the number of robots of the average total simulation time by the ratio of robots in group A.

Figure A5. Comparison of the total simulation time by the number of robots for the experiments with MT using SQF as the group B’s algorithm and TRVF and NC as the group A’s algorithm by 10% of the robots for (a) holonomic and (b) non-holonomic robots, 20% for (c) holonomic and (d) non-holonomic robots and 30% for (e) holonomic and (f) non-holonomic robots.

Figure A5. Comparison of the total simulation time by the number of robots for the experiments with MT using SQF as the group B’s algorithm and TRVF and NC as the group A’s algorithm by 10% of the robots for (a) holonomic and (b) non-holonomic robots, 20% for (c) holonomic and (d) non-holonomic robots and 30% for (e) holonomic and (f) non-holonomic robots.

Figure A6. Comparison of the total simulation time by the number of robots for the experiments with MT using SQF as the group B’s algorithm and TRVF and NC as the group A’s algorithm by 40% of the robots for (a) holonomic and (b) non-holonomic robots, 50% for (c) holonomic and (d) non-holonomic robots and 60% for (e) holonomic and (f) non-holonomic robots.

Figure A6. Comparison of the total simulation time by the number of robots for the experiments with MT using SQF as the group B’s algorithm and TRVF and NC as the group A’s algorithm by 40% of the robots for (a) holonomic and (b) non-holonomic robots, 50% for (c) holonomic and (d) non-holonomic robots and 60% for (e) holonomic and (f) non-holonomic robots.

Figure A7. Comparison of the total simulation time by the number of robots for the experiments with MT using SQF as the group B’s algorithm and TRVF and NC as the group A’s algorithm by 70% of the robots for (a) holonomic and (b) non-holonomic robots, 80% for (c) holonomic and (d) non-holonomic robots and 90% for (e) holonomic and (f) non-holonomic robots.

Figure A7. Comparison of the total simulation time by the number of robots for the experiments with MT using SQF as the group B’s algorithm and TRVF and NC as the group A’s algorithm by 70% of the robots for (a) holonomic and (b) non-holonomic robots, 80% for (c) holonomic and (d) non-holonomic robots and 90% for (e) holonomic and (f) non-holonomic robots.

Figure A8. Comparison of the total simulation time sum over the number of robots by the percentage of robots in group A using NC and TRVF as group A’s algorithms. The swarm control algorithm of the robots in group B is SQF. (a) Holonomic robots. (b) Non-holonomic robots.

Figure A8. Comparison of the total simulation time sum over the number of robots by the percentage of robots in group A using NC and TRVF as group A’s algorithms. The swarm control algorithm of the robots in group B is SQF. (a) Holonomic robots. (b) Non-holonomic robots.

From Figure A7, for 90%, TRVF is significantly better only from 200 to 300 in the holonomic case and from 180 to 300 robots in the non-holonomic cases, respectively. When using the TRVF as an group A algorithm, the lanes near the leaving trajectory of the SQF algorithm create congestion. However, with NC, they tend to go directly to the target, opening space near the SQF leaving area.

In addition, from Figure A8, NC as the group A algorithm takes less or equal time than using TRVF on average for percentages from 10% to 80% for both types of robots with statistical significance.

Appendix A.3. Selection of Alg = NC

Based on the results of Appendix A.1 and Appendix A.2, although NC is the slowest algorithm for handling the common target problem for numerous robots [8], it was the best group A algorithm in most cases with statistical significance. Therefore, NC was selected as the default group A algorithm in the experiments. In addition, MT are expected to be used for small ratios of robots in group A when assuming it would be made of auxiliary robots not knowing the control algorithm of the other group in an ad hoc teamwork scenario. Remember, however, that SQF and TRVF still outperform NC for this problem as the swarm control algorithm [23], and the results only relate to the choice of group A’s algorithm. Nevertheless, if the proportion of robots in group A and an estimated number of robots are known beforehand, the group A algorithm can be chosen based on the graphs.