Scale-Free Features in Collective Robot Foraging

: In many complex systems observed in nature, properties such as scalability, adaptivity, or rapid information exchange are often accompanied by the presence of features that are scale-free, i.e., that have no characteristic scale. Following this observation, we investigate the existence of scale-free features in artiﬁcial collective systems using simulated robot swarms. We implement a large-scale swarm performing the complex task of collective foraging, and demonstrate that several space and time features of the simulated swarm—such as number of communication links or time spent in resting state—spontaneously approach the scale-free property with moderate to strong statistical plausibility. Furthermore, we report strong correlations between the latter observation and swarm performance in terms of the number of retrieved items.


Introduction
Advances in computation have made it possible to record, simulate, and analyze multi-agent complex systems in nature, such as fish schools, bird flocks, locust swarms, and ant colonies.In many of these collective systems, various attributes were found to be scale-free [1], i.e., the attributes do not have a characteristic size or value.Examples of such scale-free features found in biological systems include, among others, (i) asymptotically scale-free correlation lengths of starling flocks [2,3]-the term asymptotic refers to the behavior of a variable (in this case spatial correlation) close to a limit (in this case an infinite flock size); (ii) scale-free fluctuations of velocity and orientation correlations in moving bacterial colonies [4]; (iii) time intervals between communication calls that follow a power law-which is the mathematical representation of a scale-free property-in pairs of zebra finches; and (iv) scale-free movement patterns found in models of foraging primates [5] or midge swarms [6].
One of the most prominent findings is that the number of interactions appear to be scale-free in various real-world networks of biological and social systems [7][8][9].Multi-agent systems benefit from scale-free communication because it enables scalable, fast and efficient information transfer [10][11][12].An essential aspect of scale-free networks is that they represent complex topologies in which only a few nodes (called hubs) have a comparably high connectivity degree [7].This small percentage of highly connected hubs makes scale-free topologies vulnerable to targeted attacks but exceptionally robust to random failures (which are likely to affect the vast majority of nodes that are not hubs) [13].Furthermore, due to the high connectivity, the network diameter is small, which means that on average, any two nodes can share their information only over a few hops [11], resulting in fast information transfer.
Inspired by the high prevalence of scale-free features in (socio-) biological systems, the aim of the current study is to examine whether scale-free attributes may also spontaneously emerge in artificial collective systems.One particularly prominent example of these systems inspired by nature is swarm robotics, where robot collaboration is an essential prerequisite for the successful execution of tasks [14][15][16][17][18][19][20].Although accurate understanding and systematic design of swarm robotic systems are considered to be among the greatest challenges of contemporary robotics [21], swarm robotics benefits strongly from the progress made in wireless communication technologies, system integration, machine learning and artificial intelligence (AI).Consequently, artificial cooperative multi-agent systems gain in importance not only in applied science and engineering but also in fundamental research, allowing the shrinkage of the reality gap of detailed modeling and the accurate simulation of distributed biological systems.
Many natural collective behaviors were modeled using artificial collective systems, including bird flocking [22], locust marching [23] or cockroach aggregation [24].However, among the most prominent challenges is the collective foraging task [18].Extensive efforts were dedicated to the study of the foraging task because of its remarkable prevalence.The foraging behavior is found in various species and subspecies across the world [25,26].In most cases, its efficient implementation is essential to group survival.
In essence, a multi-agent system performing the foraging task has the goal of collectively retrieving information and other resources from the environment.Different to a single-agent implementation, the benefit of a multi-agent system is its ability to share the accessed resources and information to enhance the overall performance.However, the multi-agent foraging task exhibits a considerable degree of complexity, making its modeling and analysis very demanding [18].Successful collective foraging often requires a delicate combination of several extensively studied multi-agent sub-behaviors such as deployment [27,28], exploration [29,30], aggregation [31,32] or information sharing [33,34].Hence, even though collective foraging itself can be considered to be a specific task within a large class of multi-agent problems, it rightfully receives separate attention in numerous contemporary studies [16,[35][36][37][38].
Moreover, collective foraging is a promising behavior for many real-world applications such as exploration by aerial vehicles [39], underwater monitoring [40], or optimization of electrical networks [41].Therefore, the foraging performance of artificial multi-agent systems, potentially in combination with other types of AI, is worthwhile investigating in depth.In particular, in robot swarms, various fundamental questions have already been addressed such as the influence of interference [42,43], regulation of information flow [33] or achievement of consensus [44].Nevertheless, other relevant questions are still open to research, including how does the distribution of individual features change in relation to the input from the environment or social interactions?Is there a connection between particular feature distributions and the performance of the swarm?To address these questions with respect to scale-free properties, we simulate the foraging behavior in a robot swarm and analyze the emergence of scale-free features.For this purpose, the complexity of the foraging task is advantageous as it offers a wide range of features that can be examined for their statistical tendency to be scale-free.
Our goal can be split into the following: (i) investigating the existence of scale-free features in a robot swarm performing the foraging task, (ii) studying the correlation between these features and the swarm performance, (iii) discussing the potential role of feedback mechanisms in the emergence of such scale-free features.
We begin with defining the robot (microscopic) and the swarm (macroscopic) behaviors in Sections 2.1 and 2.2, respectively.The link between these two levels of behaviors is formulated using statistical distributions and elaborated on in Section 2.4.In Section 2.5, we describe the experimental setup.Thereafter, in Section 3 we demonstrate the occurrence of scale-invariant features-such as those related to the communication degree or times spent in foraging or resting-and their correlation with swarm performance.Furthermore, we discuss the present feedback mechanisms that may support the emergence of scale-free features in a sophisticated set of scenario configurations.Lastly, the paper is concluded in Section 4.

Robot Behavior
We focus on the robot's decision-making process that is defined by the robot's interactions and the robot's individual preferences.The robots are situated in an arena which consists of a nest and a foraging area.Each robot can switch between two states: resting and exploring.In biology, similar behavior called forager activation has been observed in harvester ants, Pogonomyrmex barbatus, and reported in several publications [45][46][47].In the exploring state, the robot moves around searching for items-which are located only in the foraging area-to retrieve to the nest.The robot moves on a straight line until it encounters another robot or a wall, in which case a collision avoidance maneuver is initiated.In the resting state, the robot rests inside the nest and only in this state it is allowed to communicate with the neighbors within its line of sight.Specifically, each robot can broadcast a message about the success or failure of its latest exploration attempt or listen to its neighbors.Received information may either increase or decrease the robot's probability to switch from resting to exploring or vice versa.Probabilities are updated continuously using fixed probability jumps-we refer to those by the term cues, as in [14].In the following, we introduce the different probability cues used in implementing the foraging behavior, in addition to the probabilities determining the switch between the two robot states.We also consider two distinct communication modes defining the duration of information exchange.

State Switching Probabilities
Following [14], there is a minimum duration θ for the robot to stay in a certain state.The purpose of having such a threshold is to ensure that robots can perform the sub-tasks associated with this state for a certain amount of time so that necessary dynamics can take place.For instance, a minimum exploring time θ e needs to be at least as long as it takes for a robot to reach the most remote items (taking into account the constant linear speed of that robot) [16].
With this in mind, let us formulate the individual response to social and environmental cues in terms of switching probabilities.We denote {i e , i r } ∈ R + 0 , R + 0 and {s e , s r } ∈ R + 0 , R + 0 as the robot's internal (i) and social (s) cues to switch to exploring (e) or resting (r) state, respectively.The probability of a robot to switch from the resting state to the exploring state is denoted by p r→e , whereas the probability to switch from the exploring state to the resting state is denoted by p e→r .
The probabilities are updated iteratively at every simulation time step in a discrete manner as in the following: where δ η (t) is the number of 'success' minus 'failure' messages received by the robot from its neighbors at every time step spent in the resting state.Additionally, the robot's own experience is characterized using δ φ (t).This is defined for every exploration attempt as follows: with t se as the time at which the robot started its current exploration.Moreover, δ φ (t) = 0 while the robot is resting.In case p e→r < 0 it is truncated to p e→r = 0 and when p e→r > 1 it is truncated to p e→r = 1; same holds for p r→e .Please note that there is a strict difference between δ η (t) and δ φ (t): δ η (t) may be non-zero only when the robot is resting inside the nest because it is computed based on the information broadcast by the neighbors which can only be received in the nest (i.e., when the robot is in the resting state).Whereas, δ φ (t) may be non-zero only when the robot is exploring-it is computed based on the robot's own foraging experience.Table 1 lists the parameters relevant for the computation of the switching probabilities.We focus on the local communication between the robots and their influence on the global swarm behavior.A common approach is to restrict robot communication only to the area within the nest.This approach is inspired by several natural systems, in which the communication takes place mainly inside the nest or hive, such as in the case of ants or honey bees [14,[47][48][49][50].Moreover, this approach accommodates two relevant properties of foraging systems: (i) it is common that the foraging area is significantly larger than the nest area, and hence, individual encountering rates outside the nest are negligibly low; (ii) high density of individuals within the nest leads to more accurate information about the environment due to the high encounter rate of individuals that explored different, distant parts of the foraging area.
Regarding particular communication strategies, it is common to let robots broadcast the last exploration result only once, namely when the robots switch to the resting state.Henceforth, we will refer to this approach as the discontinuous communication mode (DCM), because after broadcasting the message once, the active communication of the robot is interrupted and is limited to listening.In contrast, we use the term continuous communication mode (CCM) to refer to the mode in which robots continue broadcasting the result of their last foraging attempt at every time step until they switch back to the exploring state.As we will see later, the difference between these two modes does not have substantial impact on swarm performance.However, it has a significant impact on the statistical distribution of various system features for which we study the scale-free property.

Swarm Behavior
At the macroscopic level, global behavior emerges as a result of complex interactions between the robots as well as between robots and their environment.The quality of such global behavior is evaluated with respect to quantifiable objectives.In the present study, we define the swarm performance in terms of three quantitative measures: 1.
the total number of collected items at the end of the experiment N coll 2.
the average number of collected items per time spent in collision avoidance ω ca = N coll T ca , where T ca is the total duration of all collision avoidance events in the swarm throughout the experiment.One collision avoidance event includes slowing down and turning around until the angle to the closest possible obstacle is > 25 • .In essence, ω ca reflects the trade-off between coordination and interference.

3.
the average number of collected items per time spent exploring ω e = N coll T e , where T e is the aggregate time that all individuals spent in the exploring state.

Measured Features
There is a large variety of features that could potentially have scale-free character in collective foraging.We investigate such features categorizing them into, space and time features.An overview of the measured features is given in Table 2. Space features are mostly related to the inter-robot communication according to their distribution in the arena.This is a distribution that changes over time while the robots are in motion.Among the space features, the robot's communication degree d is the most important and evident.It is defined as the number of communication links to neighbors within the robot's communication range.However, in dynamic topologies-where robots move around and neighbor lists are constantly updated-the communication degree d changes frequently.Hence, we track additionally the change of the communication degree ∆d of a robot whenever fellow robots enter or leave its communication range.Beside the communication degree, we analyze space features that reflect the foraging progress such as the difference between the number of received success and failure messages, denoted by the critical degree d c,rec .Similarly, we include features that reflect the success-degree of a particular individual by measuring the difference of success to failure messages sent by that individual, denoted by d c,sent .With respect to time features, we note that in swarm robotics the individuals are commonly subject to physical interference.Robots interfere with each other or with obstacles as a result of finite-size effects influencing the dynamics of the collective behavior [42,43,51].Therefore, we investigate the time spent on collision avoidance, denoted by τ ca .Additionally, we study time features that are related to the robot's exploring time τ e .This time can be split into foraging time τ f , i.e., the time spent on searching for items, and homing time τ h , i.e., the time spent on returning to the nest.While a long foraging time effectively increases the probability of finding items, long homing times indicate overcrowding close to the nest.Finally, another relevant time feature is the resting time τ r that includes the duration of robot interaction within the nest.

Data Analysis
In complex systems such as swarm robotics, the statistical analysis of relevant system properties paves the way to mathematical modeling, useful simplifications, or inference of long-term behaviors.Consequently, it helps in defining the link between the individual robot behavior and the emergent global swarm behavior, referred to as the micro-macro link [52].In our study, we focus on how the collective foraging behavior can be related to the scale-freeness of a set of individual and global features.The main statistical characteristic of scale-free features is that they are distributed according to a power law [1,53].Therefore, to identify scale-free features in our simulated swarms, it is of central importance to measure the statistical distribution of these features and to perform a sound power law fitting procedure.

Power Law Fitting Procedures
To verify whether a feature is scale-free, we use a set of techniques that are described in [53][54][55] for fitting its distribution by the power law distribution.The power law distribution takes the form of a straight line on a log-log scale of p(x).However, most real-world data displays significant fluctuations due to randomness.When fitting power law to the data, random fluctuations are considered by the statistical value p which represents the goodness-of-fit.When p < 0.1 the power law fit can be considered to be unreliable [54].Furthermore, power law behavior emerges mostly only in the tail of the distribution, i.e., for higher values of x above a statistically determined lower bound x min [53].Please note that this effectively reduces data set to fit by the power law, which is important to keep in mind by considering a ratio of the total number of data points to the points that satisfy the condition x > x min .Finally, there are several other statistical distributions that may resemble the characteristic straight-line tendency of a power law on a log-log plot.Hence, for a sound statistical analysis it is important to compare the power law fit to other statistical models [54][55][56].More precisely, the power law fitting procedure can be summarized by the following three steps: 1.
Using maximum likelihood estimation, fit the data by the power law distribution where α is the scaling parameter and x min is the lower bound.In particular, α and x min are estimated using procedures described in [54].

2.
Apply Kolmogorov-Smirnov statistic to carry out the goodness-of-fit tests and verify the above results.
Here, essentially, a large set of synthetic data is generated from a power law distribution (with α and x min found in 1.) and their distances to their respective fits are compared to the distance between the empiric data and the best-fit found in step 1.The outcome of this procedure is the p-ratio which estimates the contribution of random fluctuations.If p < 0.1 the power law fit found in step 1 is very likely to be due to inherent randomness in the empiric data.Moreover, the p-value is unreliable when the data set is too small.Therefore, the percentage of data points N data,pl that lie above x min should be 10% or higher.

3.
Finally, even if p > 0.1, the power law fit might be not the only model that fits the data well.Consequently, complete the above steps for a set of other potential distributions including exponential or lognormal distributions.Then, compare the resulting best fits to the one obtained for the power law by computing the ratio R. The latter is defined as the log likelihood of the power law over the log likelihood of another distribution.If R > 0, power law is the statistically superior fit.Although this last step still does not give us the certainty that the data is power law distributed, it makes the hypothesis more plausible.For this step we used an open-access Python toolbox [56].

Quality Ratio ρ q
Given a high quantity of empiric data sets, it is useful to find an automated way for the evaluation of the power law fits.For the analysis of our experiments, we introduce a quality ratio ρ q which we use as a practical estimate of the plausibility of a (truncated) power law fit based on the well-known rigorous statistical tests described above.The quality ratio ρ q includes the three criteria discussed in Section 2.4.1:p-value, N data,pl and the number of likelihood-ratio tests resulting in R > 0. We account for these criteria by defining ρ q as the product of ρ q,p , ρ q,data and ρ q,lrt : 1.
First, we begin with the p-value.As mentioned above, the linear shape of the data distribution on a log-log plot can be mainly attributed to random fluctuations if p < 0.1.Taking this into account, we design ρ q to be a binary piecewise function evaluating the goodness-of-fit in terms of the p-value: This way, we take into account the possibility that random fluctuations may be present but as soon as p > 0.1 we do not assign the precise value of p to the ranking of the fit.The reason is that random fluctuations might be present even if the data is in fact power law distributed.In that case the p-value could be very low even if the data is in fact power law distributed.In general, it might be more substantial to consider the size of the fitted data set and to compare the power law fit to other important distributions [54,56].2.
Second, ρ q,data , denotes the ratio of the data which is fit by the (truncated) power law N data,pl to the total number of data points N data,tot .
Third, ρ q,lrt represents the fraction of likelihood-ratio-tests in which the (truncated) power law fit proved to be statistically more plausible than other distributions.To include the quality of the (truncated) power law fit as compared to other distributions, we count how many times n lrt,pl we obtained R > 0 from the likelihood-ratio tests.We compare the power law fit to six distributions: truncated power law, exponential, stretched exponential, lognormal, positive lognormal and normal; all of them are implemented in [56] (except the normal distribution).Hence, we use the piecewise function, , otherwise.
where we added 1 to n lrt,pl to account for the possibility that the likelihood-ratio test yields R ≈ 0, in which case the support for the power law fit is neither strengthened nor weakened.Please note that in Equation ( 7) we set ρ q,lrt = 0 if at least one distribution is a more reliable model than the power law.However, it is important to remember that our simulated systems are meant to include real-world attributes (e.g., finite-size effects, physical interference, line-of-sight interruptions during communication) and therefore deviate from ideal systems.Consequently, the assumption of power law (i.e., scale-free) distribution might be distorted and needs to be corrected.The deviation is often particularly distinct in the heavy tail.Therefore, one common correction technique is to consider the power law distribution with an exponential cutoff (also known as truncated power law) [57]: where λ is the scaling parameter of the exponential decay and Γ(y, z) is the upper incomplete gamma function.While Equation ( 4) directly implies that the feature is scale-free, Equation ( 8) describes an asymptotic scale-freeness in the limit λx → 0. This equation approaches the power law distribution asymptotically for λx → 0 and the exponential distribution for xλ 1, respectively.Thus, accepting that our systems are significantly constrained within physical boundaries we can slightly soften the criteria given by Equation (7) in the following way: If the truncated power law passes more likelihood-ratio tests than the power law fit, i.e., if n lrt,tpl > n lrt,pl , we consider the success-ratio of the former.In short: Finally, including all the above criteria, we define the quality ratio: Consequently, we obtain ρ q = 0 if p ≤ 0.1 or R <0.Conversely, ρ q = 1 in the case of p > 0.1, N data,pl = N data,tot and n lrt = 6, which is an unlikely but nevertheless possible scenario.Using this ranking, we can link the quality of a fit to a quantifiable value and describe the support for the (truncated) power law as illustrated in Table 3.
Table 3. Classification of the power law fit quality with respect to the quality ratio ρ q used in our study.

< ρ q
The denominator value represents the total number of considered distributions (i.e., the power law and the six alternative distributions we compare it to).The lower limit of the 'moderate' classification corresponds to the case with ρ q,p = 1, ρ q,data = 0.1 and ρ q,lrt = 1 7 -i.e., at least 10% of the data is included in the fit and none of the alternative distributions is a statistically better fit than the power law.The upper limit considers the case with ρ q,p = 1 and ρ q,data • ρ q,lrt = 0.5 7 -i.e., either the fit includes a high number of data or the power law is statistically a better fit than other distributions.Please note that ρ q multiplicatively combines standard power law fitting techniques [53][54][55][56] into a quantitative estimate of the quality of the (truncated) power law fit.
It is important to emphasize that even if the hypothesis of the data following the power law distribution is found to be plausible using the above statistical analysis, care needs to be taken when interpreting this observation.Firstly, there is still no guarantee that the data is in fact power law distributed and although our rigorous analysis includes several common distributions, other non-obvious distributions may prove to be a better fit.Secondly, the power law fit may be valid only for a small fraction of data.However, as the power law behavior is commonly found for a subset of data, namely at the tail of the distribution, the group that displays power law (i.e., scale-free) behavior includes individuals that stand out from the rest of the swarm by having features with values that are significantly above average.The way in which such individuals impact the global swarm performance remains an open question worth investigating.

Correlation Measures
To examine the presence of correlations between the support for the power law distribution (i.e., the value of ρ q ) and the swarm performance it is important to use an appropriate correlation measure.One of the most prominent correlation measures is the Pearson correlation coefficient [58,59].It evaluates the quality of a linear association between two distributions.In essence, it calculates the covariance of the mean values of two distributions, over the root of their standard deviations.It is closely related to linear regression and does not require the data to be normally distributed.Despite its mathematical simplicity it is an appropriate correlation measure for many distributions and, therefore, is widely used [60][61][62][63].
However, one could argue that the Pearson correlation coefficient is not ideal for skewed distributions with strong outliers.Popular alternatives are the Spearman's rank and the Kendall's tau correlation coefficients [62][63][64][65].Both are based on generating ranked distributions by assigning a rank to each variable with respect to its value.The correlation coefficient is then given as a measure of the association between the two ranked distributions.Consequently, both correlation metrics are robust to outliers and suitable for non-linear distributions.
Although both correlation measures commonly return very similar results, Kendall's tau handles ties (i.e., cases in which there is no difference between the ranks) in a mathematically more straightforward way.More precisely, Kendall's tau returns the density difference between concordant and discordant pairs.Consider two vectors of length n, (x 1 , x 2 , ..., x n ) and (y 1 , y 2 , ..., y n ).Concordant pairs are pairs of data points that satisfy sgn(x i − x j )sgn(y i − y j ) > 0 (where sgn (z) is a sign-function equal to +1 if z > 0, −1 if z < 0 and 0 if z = 0); similarly, discordant pairs satisfy sgn(x i − x j )sgn(y i − y j ) < 0. Furthermore, ties are pairs for which x i = x j or y i = y j .Hence, with n c (n d ) as the number of concordant (discordant) pairs, respectively, and n x (n y ) as the number of ties in x (y), respectively, the Kendall's tau (also known as the Kendall's tau-b) is given by [66]: with where and similarly, for δ i,j and δ (y) i,j .

Simulation Setup
We designed and implemented a set of physics-based simulations using the state-of-the-art simulator for large-scale swarms, ARGoS [15].An overview of all parameter values used in our simulations is given in Table 4.The simulations are conducted in a square-shaped arena, which is confined within four walls, each being of the length of 50 m.The arena is divided into two regions: (i) the nest A n : it is the gray 10 × 50 m 2 area in Figure 1a, and (ii) the foraging area A f : it is the white 40 × 50 m 2 area in Figure 1a.
The items are scattered uniformly over the foraging area, and keep reappearing after robots retrieve them to the nest-as in [14]-with constant probability.This prevents the system from drifting into an absorbing state in which there are no items left to recover.A phototaxis behavior is used to assist the robots in leaving and re-visiting the nest.For that purpose, light beacons are positioned equidistantly at the nest wall (yellow dots at the bottom of Figure 1a).Their light is perceived by the robots' light sensors.Each robot is programmed to move away from the beacons when it needs to leave the nest and towards the beacons when it needs to return.We use a homogeneous swarm of Footbots (see http://www.swarmanoid.org/swarmanoid_hardware.php) in our simulations, and the communication radius of the robots is set such that the fraction of the circular communication area around the robot is 0.982 % of the nest area, which is close to the fraction used in [14].For better readability, we will limit our discussion of the robot states to only resting and exploring.While the former is distinct, the latter is composed of further states of which only foraging and homing are relevant because they are the most time consuming (for a detailed list of the robot states please see Supplementary Material Section 2).
At the beginning of the simulation, each robot switches from resting to exploring with a probability of 0.01.Consequently, within the first 500 time steps (ts) most robots leave the nest.After another ≈ 500 ts most of the swarm returns, with or without an item.Although this behavior is subsequently repeated several times, the number of simultaneously switching robots gradually decreases, and the switching rate from resting to exploring (or vice versa) approaches a constant limit.In most cases, the system approached an equilibrium after 5 • 10 3 ts (for the given arena, item density, and swarm size).Our measurements of the system features begin from that time instance on-wards and the experiment proceeds for another T = 10 4 ts.
Furthermore, to conduct a solid statistical study we use large-scale swarms with N robots = 950 units which is up to an order of magnitude higher than what is commonly used [14,16,36,42].We selected the value of N robots by running preliminary experiments, in which we observed for this particular swarm size-under the given arena and item density-a maximum in swarm performance.
Finally, in our experiments, the most important means of influencing the swarm dynamics is by adjusting the numerical values of the internal and social cues-i e , i r and s e , s r , respectively-at the start of each experiment.We consider a spectrum of 256 distinct scenario configurations, which differ by the 4-tuples a = (s e , i e , i r , s r ) drawn from: Ω := {a : s e , i e , i r , s r ∈ {0.0, 0.01, 0.5, 0.9}}.
The rationale behind the choice of these parameter values is to include four fundamentally different kinds of cue impact on swarm dynamics: (i) none (ii) low (iii) intermediate and (iv) high.Please note that any additional value in the set a greatly increases the associated computational and analytic effort-as the number of scenarios scales with dim(a) 4 .However, based on preliminary results, additional values would offer potentially little informative gain (at the current stage) because the swarm dynamics would be similar to a mix of the dynamics generated by the above values.

Results and Discussion
We performed simulations with all combinations of cues and communication modes.Each simulation was repeated with 30 random seeds and the data analysis procedure was carried out as discussed in the previous section.

Presence of Power Law Distributed Features
The analysis of our simulation data shows that in most scenarios there was only weak or no statistical support for the (truncated) power law distribution (see Figure 2A).In particular, in roughly half of all scenarios no power distributed features were found.This observation suggests that in the present system, scale-free features are rare.Nevertheless, we found 245 + 71 = 316 (truncated) power law distributions with moderate or strong statistical plausibility for different features in various scenario configurations and for both communication modes, DCM and CCM.Thus, our findings are in line with a recent study showing that scale-free networks may occur rarely but across different areas [55].
As the scatter plots in Figure 2 show, most of the distributions with weak or moderate support for power law are concentrated below or close to the average values of swarm performance while the distributions with strong support for power law are associated with above-average performance in terms of N coll and ω ca .Swarm performance was measured using (i) the number of items retrieved by the robots, N coll , (ii) the average number of collected items per time spent on collision avoidance ω ca , and (iii) the average number of collected items per time spent exploring ω e .Figure 3 shows the values recorded for these three metrics under both continuous and DCMs and for the entire range of 256 scenario (cues) configurations, respectively.Repeating performance patterns can be observed over different sets of configurations.The regions over which these patterns emerge are (from left to right): (i) all scenarios with s e = 0 and i e = 0, i.e., constant p r→e (blue region with a left tilted mesh in Figure 3), (ii) all scenarios with s e = 0 and i e > 0, i.e., no social and only internal influence on p r→e .This region is henceforth referred to as NS e (shown in orange, no mesh), (iii) all scenarios with s e = 0.01, i.e., low social impact on p r→e (green region with vertical mesh, henceforth denoted as LS e ), and (iv) all scenarios with s e = 0.5 or s e = 0.9, i.e., high social influence on p r→e (red region with right tilted mesh, henceforth denoted as HS e ).Please note that in all four regions p e→r is altered in the same way, i.e., for s r and i r all values from {0.0, 0.01, 0.5, 0.9} are included.The best swarm performance in terms of N coll and ω ca emerges when the influence of internal cues on the swarm dynamics is negligible compared to social cues, i.e., when  3. (b-d) Log-linear scatter plots relating the power law fit quality ratio ρ q to the swarm performance in terms of N coll , ω ca and ω e , respectively.The vertical dashed lines indicate the mean performance values while the horizontal dashed lines separate the quality categorizations taken from Table 3.  c,d) ω ca and (e,f) ω e , respectively.For each performance measure, 256 scenario configurations were implemented (i.e., with all cue values from Equation ( 14)), using one of two communication modes: DCM (left) and CCM (right).The x-axis represents the IDs of the scenario configurations.The colors and the mesh patterns highlight regions that display different dynamics.Apart from (f), in all plots the red dots mark the scenarios in which the feature mentioned in the inset demonstrated a high value of ρ q , i.e., there was a strong support for the distribution to be power law.In (f), the red dots mark the scenarios with moderate support.See Supplementary Material Section 3 for combined plots of N coll and d distribution in CCM over the complete set of 256 scenario configurations.
The best performance levels in terms of N coll and ω ca were reached over the LS e region.For instance, the maxima of N coll and ω ca correspond to the scenario configurations in which s e = 0.01, i r = 0 and s r ≥ 0.5.For the same configurations, (truncated) power law distributions of space features were found in the CCM (examples shown in Figure 4).Contrary to CCM, in the DCM the robot interactions are interrupted.These interruptions may explain why, in DCM, space features such as communication degree tend to not follow a power law distribution (weak overall support for the presence of a power law behavior).Nevertheless, we found fits with moderate to strong support for (truncated) power law to time features, such as τ r and τ ca , demonstrated in Figure 5.The best power law fits of the DCM correspond to the peaks in swarm performance in terms of N coll and ω ca over the HS e regions.c,d).These scenarios belong to a subset of the best swarm performance with respect to ω ca .Please note that λ = 0 for the fits of τ ca , indicating better support for the power law fit than for truncated power law.
The third performance measure, i.e., ω e , reached its best values over the NS e region.Its maxima correspond to cases where s e = 0 and s r = 0. Interestingly, for these scenario configurations we found fits with moderate to strong support for (truncated) power law to the data of ∆d, i.e., the change of the average communication degree of the robot (examples shown in Figure 6).This is an interesting finding because it indicates that a communication feature may be power law distributed also in those scenarios in which the swarm tries to minimize the number of foraging robots and maximize the number of resting ones.Moreover, in most ∆d distributions with strong or moderate support for the power law, the fit includes only 10-20% of data points.The reason for the relatively low ratio of power law fitted data is that the tail of the distribution is likely to represent by the fraction of robots that rest or move close to the border between the nest and the foraging area.
In general, the findings suggest that internal cues (in the absence of social cues) keep robots at the edge of minimal activity while social cues (in the absence of internal cues) drive the robots towards maximal activity.

Correlation with Swarm Performance
In the previous section we have illustrated that swarm performance is likely to reach its peaks over cue configurations that include asymptotically scale-free space or time features (see Figure 3).In this section, we analyze this observation statistically, using correlation measures such as the Pearson and Kendall's tau rank correlation coefficients introduced in Section 2.4.3.However, note that both correlation measures have strengths and shortcomings.On the one hand, while the Pearson correlation coefficient is widely used and has an elegant mathematical form, it is sensible to outliers and may not be appropriate for non-linear distributions.On the other hand, Kendall's tau is suitable for non-linear distributions as well as robust to outliers.Nevertheless, reducing the values to ranks may disregard the significance of the variable's value being far from the average.In particular, replacing the real value of the quality ratio ρ q by its rank leads to loss of information about the extent to which ρ q represents the quality of the power law distribution.Moreover, following the definition of Kendall's tau in Equations ( 11)- (15), each difference between data point pairs is assigned the same weight which may not always be appropriate.For instance, consider the ranked swarm performance in terms of N coll and the corresponding distribution of ρ q in CCM for d in Figure 7a and for d c,sent in Figure 7b, respectively.In both cases, Kendall's tau defined by Equations ( 11)-( 15) returns values indicating no correlation (i.e., τ Kendall = 0.02 and τ Kendall = −0.04,respectively).However, as evident in Figure 7, both cases show different dynamics, with ρ q for d following N coll more closely than for d c,sent .The main reason is that the dominant fluctuations of ρ q close to zero are assigned the same weight (i.e., rank step 1) as the more permanent increase of ρ q for high values of N coll .Similar considerations hold for the other features and the DCM.To account for this type of behavior, we use a generalization of Equation ( 15) that weights the ranking steps by a parameter κ, which is relative to the average change, such that: A high correlation coefficient means: the better the quality of (truncated) power law distribution, the higher the likelihood that the swarm performed well.Cells highlighted in blue show moderate or strong correlations.
The correlation coefficients confirm the observation, supported by the data shown in Figures 2 and 3, that most power law distributions with strong support (i.e., high ρ q ) appear in scenarios with peak swarm performance.To further examine this observation, we consider the correlations of the swarm performance with different ρ q support classifications (based on Table 3).As Table 6 shows, there are moderate and strong positive correlations between features with strong support for power law distribution and swarm performance in terms of N coll and ω ca for both communication modes.This suggests that the observation of scale-free features is more likely in scenarios in which the agents are more successful in retrieving a high number of food items.

The Role of Feedback Mechanisms in the Emergence of Scale-Free Features
An attribute of complex systems that is widely known to support the emergence of scale-free characteristics is the presence of (positive and negative) feedback loops [1,53,71].We specify the feedback effect to be positive or negative based on the individual response to the information input from its neighborhood.Hence, we refer to the feedback mechanism as positive feedback if it pushes the individuals to the same state as the state of the majority, whereas negative feedback pushes them away from it.
Most scale-free features were found in scenarios in which (i) the robot behavior was dominated by social interactions, (ii) the swarm attempted to balance positive and negative feedback and (iii) the swarm displayed a tendency towards active exploration.In particular, in CCM, the first 17 scenarios sorted by ρ q in descending order were found over the LS e region and with i r = 0. Similarly, in DCM, the first 28 scenarios were found over the HS e region and with i r = 0. To understand this, it is necessary to consider in more detail the impact of each cue on swarm dynamics and the feedback mechanisms.
For conciseness, we focus our system analysis on CCM and its most relevant set of parameter configurations.Similar conclusions hold for DCM.In particular, we can simplify our analysis based on the repeating patterns of swarm performance (see Figure 3) and the following observations: (i) The swarm performance is qualitatively very similar between the cue values 0.5 and 0.9 (for all four cues).Thus, we focus, in the following, on {a : s e , i e , i r , s r ∈ {0.0, 0.01, 0.5}}.(ii) The cue i e has a negligible impact on the foraging dynamics when s e > 0. By neglecting scenarios in which i e = 0, except those with s e = 0, we can further shorten the set of relevant scenarios.Finally, (iii) there are significant differences in the dynamics between scenario configurations with i r = 0 and those with i r > 0 but negligible differences between i r = 0.01 and i r = 0.5, 0.9.Thus, we focus, in the following, only on scenarios with either i r = 0 or i r = 0.01. Figure 8 shows the final set of 24 scenarios relevant to the discussion below.Please note that (ii) and (iii) are consequences of the internal cues i e and i r acting only on exploring robots.In the exploring state, the crucial parameter is p e→r because it defines the probability to stop exploring and change to resting.A non-zero value of the internal cue i r has a substantial impact on dynamics as it alters p e→r after each exploration attempt.As the likelihood of finding and retrieving a food item is low, i r mostly reduces p e→r .The more p e→r is lowered by i r , the less likely the robot is to find a food item during the next exploration attempt.Thus, i r has a strong inhibitory influence on the swarm's exploration activity.Consequently, there is a significant difference in swarm performance between the scenario configurations with i r = 0 and those with i r > 0. As the swarm actively attempts to explore the environment and collect food items, the influence of i r > 0 can be considered an important driver of negative feedback.In contrast, p r→e acts only on resting robots.Consequently, any change of p r→e through i e is easily distorted by s e , i.e., the social interactions with the neighborhood of the resting robot.Hence, only when s e = 0 there is an inhibitory impact of i e on swarm dynamics (similar to i r , due to the scarcity of food items), otherwise i e is negligible.In short, when the probability of finding a food item is low, with i r = 0 and s e > 0 enabling the swarm to significantly damp the feedback mechanisms that drive the swarm towards inactivity.
Next, we consider the particular contributions of social cues s e and s r .Social interactions represent a direct form of feedback loops, enabling the swarm to drift towards an absorbing state (e.g., uninterrupted resting or exploring) or maintain a balance between positive and negative feedback.In general, note that high values of s e often lead to p r→e = 0 due to the high probability of encountering a robot with a failed exploration (due to the low density of food items).With such relatively high values of s e , p r→e can be reduced to zero within a few time steps.By contrast, with s e = 0.01, p r→e does not fluctuate as strongly.Similar considerations hold for s r and p e→r .In terms of active exploration, i.e., long exploration times and high number of retrieved items, it is beneficial for the swarm to have robots with p r→e = 1 and p e→r = 0. Indeed, in the present system we observe that the swarm approaches such behavior for s r = 0.5 and s e = 0.01 (i.e., over the LS e region).More importantly, such balance of s r and s e allows positive and negative feedback loops to coexist with the positive feedback being slightly more dominant.Due to this feedback coexistence, a robot that happened to be surrounded by unsuccessful neighbors will tend to have low p r→e and high p e→r , i.e., its resting time τ r will increase (together with its d or d c,rec ) and vice versa.Over time, such dynamics will result in robots that are increasingly inactive (with increasingly higher τ r , d or d c,rec ) and robots that are increasingly active (with increasingly lower τ r , d or d c,rec ).When the majority tends towards active exploration, the inactive group of robots experiences negative feedback and while the active group is subject to positive feedback.The prevalence of the positive feedback decreases the number of consistently resting robots significantly below the number of consistently exploring ones.Ultimately, this leads to skewed or heavy tailed distributions, such as the power law and, consequently, to the emergence of scale-free features.Similar considerations apply to the DCM over the HS e region.The difference is that in DCM each robot can broadcast its exploration result only once.Thus, s e needs to have high values for dynamics similar to CCM to emerge.
To illustrate the above considerations, let us examine the scenario configuration of CCM with s e = 0.01, i e = 0.9, i r = 0.0, and s r = 0.5 (see Figure 4a), in which a high performance value of N coll was observed (i.e., this scenario is similar to the peak in Figure 8).The value of s r has a high impact on p e→r (the probability to transition to resting).For example, if a robot receives at least two 'success' messages more than 'failure' messages-i.e., if δ η (t r ) ≥ 2 in Equation (2)-its p e→r drops to zero.When p e→r = 0 and i r = 0.0, the robot will stop exploring only if it finds an item.During the subsequent resting, this robot is likely to cause one of its neighbors to reach p e→r = 0, which repeats an analogous cycle of events.The corresponding dynamics can be translated in terms of the positive feedback pushing the robots out of the nest and increasing the number of robots in the foraging area (i.e., in the exploring state).In the long term, due to the positive feedback, the swarm drifts towards the absorbing state in which all robots have p e→r = 0 and p r→e = 1.0.In the short term, while most robots is exploring, some robots remain in the nest, e.g., due to crowding at the entrance of the nest.Those robots have a higher number of neighbors because the nest is significantly smaller than the foraging area.Therefore, during this crowding behavior, the swarm experiences the coexistence of positive and negative feedback loops.A specific balance between these feedback loops may lead to the emergence of scale-free features such as the space feature d (for which, indeed, the above mentioned scenario configuration has one of the best truncated power law fits with ρ q ≈ 0.23, shown in Figure 4a).Similar considerations hold for other CCM examples presented in Figure 4 or the DCM examples in Figure 5.
The above example demonstrates positive feedback regarding the exploring state.However, under some configurations, positive feedback can also be observed around the resting state.For example, during the crowding behavior in the nest, a robot which is surrounded by a high number of resting neighbors is likely to get 'stuck' and be unable to leave the nest.This robot will eventually switch to the resting state and broadcast a 'failure' message.Neighbors that receive this message will decrease their probability to explore p r→e , and, through physical interference, lower the neighboring robots' chances to leave the nest.Consequently, positive feedback at the resting state may lead, in the long term, to an increase of the average communication degree of the resting robots and the emergence of power law distributed features, alongside the occurrence of outstanding robots whose features such as ∆d (examples shown in Figure 6), τ r or τ ca exhibit exceptionally above-average values.

Conclusions
In this paper, we have investigated the interplay between scale-free features and swarm dynamics of a foraging swarm.Our results demonstrate that in the studied system (i) several space and time features tend to be asymptotically scale-free for multiple parameter configurations; (ii) the emergence of scale-free features can be attributed to the presence of positive/negative feedback mechanisms.Furthermore, (iii) in several cases the swarm performance is moderately or strongly correlated with the tendency of space and time features to follow the power law distribution-which is the mathematical backbone of the scale-free property.
This study serves as a first step towards a better understanding of the interplay between the presence of scale-free features and the swarm behavior in terms of collective performance.Although our results do not indicate a causal relationship, we found conclusive evidence for a close connection between scale-free features and swarm performance.Moreover, our analysis of power law distributed features shows a strong link between the microscopic behavior of robots determined by specific cues and the macroscopic behavior of the entire swarm exhibiting peak performance.However, care needs to be taken when considering cases where the power law fit includes only a small fraction of data, as focusing on a small subgroup that plausibly displays scale-free features may disregard a significant piece of information about the global swarm behavior.
Please note that the presented exploratory study was conducted with an emphasis on whether scale-free features may emerge autonomously in artificial multi-agent systems, without focusing on why they do so.Hence, more sophisticated work is needed to precisely understand the exact causes for the emergent scale-free characteristics in our systems.For instance, strong feedback mechanisms may push the system close to a critical point at which a phase-transition occurs.In case of a continuous phase-transition, the latter is known to be associated with the emergence of scale-free features [1,53].In fact, using approximations (such as the assumption of a well-mixed system) it could be shown analytically that the social or internal cues can be used as control parameters, moving the system between its phases (e.g., phases in which the number of resting robots is minimized or maximized).However, if the system approaching a phase-transition is the cause for the emergence of scale-free features in our experiments, we expect to find a correlation length (i.e., the distance over which one robot influences another, directly or indirectly) that is longer than the size of the system (e.g., the length of the nest).In contrast, our preliminary analysis indicated the opposite behavior: as we approached those scenarios in which scale-free characteristics were observed, the correlation length decreased below the system size.In general, a detailed finite-size-scaling analysis is necessary to explain our findings more thoroughly as well as reveal which (physical) boundaries are most relevant and what impact they have on the system dynamics.
The canonical foraging task continues drawing scientific attention due its importance and prevalence in nature as well as artificial systems.In addition, the complexity of collective foraging as a combination of several sub-behaviors allows the modeling and analysis of a large number of scenarios and examine various features.For these reasons, we focused exclusively on the foraging behavior.However, it would be interesting to extend the scope and investigate other multi-agent tasks, such as aggregation or flocking,

Figure 1 .
Figure 1.Illustrations of the arena.(a) A snapshot from a simulation in ARGoS.Gray area: nest; white area: foraging field; black dots: items; blue objects: Footbots; light-blue lines: communication (range-and-bearing) links.(b) 3D view on the same arena.In both figures, the communication links are formed only for resting robots inside the nest, as in our experiments moving robots neither broadcast nor listen to any messages.

Figure 2 .
Figure 2. (a) Feature data sets obtained from simulations, sorted by the type of statistical support for a corresponding power law fit.The classifications follow Table3.(b-d) Log-linear scatter plots relating the power law fit quality ratio ρ q to the swarm performance in terms of N coll , ω ca and ω e , respectively.The vertical dashed lines indicate the mean performance values while the horizontal dashed lines separate the quality categorizations taken from Table3.

Figure 3 .
Figure 3. Swarm performance in terms of (a,b) N coll; (c,d) ω ca and (e,f) ω e , respectively.For each performance measure, 256 scenario configurations were implemented (i.e., with all cue values from Equation (14)), using one of two communication modes: DCM (left) and CCM (right).The x-axis represents the IDs of the scenario configurations.The colors and the mesh patterns highlight regions that display different dynamics.Apart from (f), in all plots the red dots mark the scenarios in which the feature mentioned in the inset demonstrated a high value of ρ q , i.e., there was a strong support for the distribution to be power law.In (f), the red dots mark the scenarios with moderate support.See Supplementary Material Section 3 for combined plots of N coll and d distribution in CCM over the complete set of 256 scenario configurations.

Figure 4 .
Figure 4. Log-log scale plots of the degree d (top) and the critical degree d c,rec (bottom) distributions in CCM.The black lines represent the corresponding truncated power law fits.The insets show the fit parameters as well as the scenario configurations.The plots (a) and (b) differ by the scenario configurations shown in the insets; similarly for (c) and (d).These scenarios are among the top five swarm performances with respect to N coll .Please note that λx min is relatively small, i.e., power law is a good fit for x close to x min .See Supplementary Material Section 3 for plots of d in CCM over all 256 scenario configurations.

Figure 5 .
Figure 5. Log-log scale plots of the resting time τ r (top) and the collision avoidance time τ ca (bottom) distributions in DCM.The black lines represent the corresponding truncated power law fits.The insets show the fit parameters as well as the scenario configurations.The plots (a,b) differ by the scenario configurations shown in the insets; similarly for (c,d).These scenarios belong to a subset of the best swarm performance with respect to ω ca .Please note that λ = 0 for the fits of τ ca , indicating better support for the power law fit than for truncated power law.

Figure 6 .
Figure 6.Log-log scale plots of ∆d per 100 ts in CCM.The black lines represent the corresponding truncated power law fits.The insets show the fit parameters as well as the scenario configurations.The plots (a,b) differ by the scenario configurations shown in the insets.These scenarios are among the best swarm performances with respect to ω e .

Figure 8 .
Figure 8. Number of collected items for a selected set of 24 scenario configurations a in CCM.The data labels show the corresponding cue values of a = (s e , i e , i r , s r ).

Table 1 .
An overview of parameters defining the switching probabilities.

Table 2 .
An overview of the investigated space and time features.

Table 4 .
Robot and arena parameters used for the simulation setup.

Table 5 .
Correlation coefficients quantifying correlations of ρ q with N coll , ω ca or ω e .

Table 6 .
Correlation coefficients between ρ q and N coll , ω ca or ω e for different power law support classifications.

for pl Correlation with N coll Correlation with ω ca Correlation with ω e
Correlation of the swarm performance with different categories of power law distribution support.Cells highlighted in blue show moderate or strong correlation values.