Characterizing The Limits of Linear Modeling of Non-Linear Swarm Behaviors

We study the limits of linear modeling of swarm behavior by characterizing the inflection point beyond which linear models of swarm collective behavior break down. The problem we consider is a central place object gathering task. We design a linear model which strives to capture the underlying dynamics of object gathering in robot swarms from first principles, rather than extensively relying on post-hoc model fitting. We evaluate our model with swarms of up to 8,000 robots in simulation, demonstrating that it accurately captures underlying swarm behavioral dynamics when the swarm can be approximated using the mean-field model, and when it cannot, and finite-size effects are present. We further apply our model to swarms exhibiting non-linear behaviors, and show that it still provides accurate predictions in some scenarios, thereby establishing better practical limits on linear modeling of swarm behaviors.


Introduction
Swarm Robotics (SR) is the study of the coordination of large numbers of simple robots (Ş ahin, 2005). SR systems can be homogeneous, i.e., single robot type and identical control software, or heterogeneous, i.e., multiple robot types and control software (Prorok et al, 2017;Ramachandran et al, 2020). The main differentiating factors between SR systems and multi-agent robotics systems stem from the mechanisms on which SR systems are based. Historically, these were principles of biological mimicry or problem solving techniques inspired by natural systems of agents such as bees, ants, and termites (Labella and Dorigo, 2006), though modern SR systems typically incorporate more mathematically rigorous elements of conventional multi-robot system design (Castello et al, 2016).
In this work, we study the central place foraging problem in which robots gather objects (blocks) across a finite operating arena and bring them to a central location (nest) under various conditions and constraints. Foraging is one of the most studied applications of SR, due to its straightforward mapping to real-world applications (Hecker and Moses, 2015); for an extensive discussion of the state of the art, see (Lu et al, 2020). The complexity of the foraging task frequently gives rise to non-linear behavioral dynamics caused by inter-robot interactions.
It has been established that above a certain swarm density ρ S (ratio of swarm size to arena size) (Sugawara and Sano, 1997;Hamann, 2013), interactions between robots can produce non-linear behaviors which are not predictable from the individual components, e.g., from the swarm control algorithm (Cotsaftis, 2009;Georgé and Gleizes, 2005;Hunt, 2020;De Wolf and Holvoet, 2005). Hence, the collective performance of a swarm S of N robots each running an identical control algorithm κ can be a non-linear function of the behavior of a system of N independent robots (Harwell et al, 2020). Below this density threshold, swarm behavior can be well approximated using linear models.
Many real problems across scales can be tackled using swarms in low density environments independently of swarm size, in which each robot is responsible for a large area on the order of 100m 2 ; these include indoor warehouses (≈ 64m 2 ), outdoor search and rescue, precision agriculture, and field monitoring (600,000m 2 ≈ 150 acres). Recent results summarizing the challenges of moving swarms into the real world argue that directing research towards low density swarm applications is critical (Tarapore et al, 2020). Some researchers hold that swarms with low ρ S are not properly characterized as swarms, and are instead systems of independent robots because they lack the high level of inter-agent interaction which characterizes natural swarms. We argue that a high level of inter-agent interaction, and therefore potentially emergent self-organization, is only one of the defining properties of swarms. Low density swarms exhibiting the properties of scalability, flexibility, and robustness can properly still be considered swarms (Harwell et al, 2020). The exact value of ρ S at which a given linear model of swarm behavior breaks down is influenced by many factors, including the control algorithm κ, the number of robots N, and characteristics of the problem being solved, so in general it cannot be determined a priori.
In this work we seek to characterize, as a function of ρ S , the practical limits of using linear models to model non-linear behaviors of a swarm engaged in real-world foraging. Specifically, we will show that non-linear behavioral dynamics in the foraging problem can be captured using linear modeling which is derived directly from problem features, such as swarm density, swarm size, and robot control algorithm.

Motivation and related work
Many SR systems have been designed around imitating natural systems, and use heuristic decision making (Castello et al, 2016) rather than combining natural principles with a strong mathematical grounding (Talamali et al, 2020). Nevertheless, heuristic approaches to swarm control have been effective for robots that operate with incomplete information and limited computing power. SR researchers taking this approach average large numbers of simulation runs to develop more accurate models of swarm behavior and obtain empirical insights into real-world problems (Harwell and Gini, 2019). The emphasis on empirical rather than mathematical models has been a chief impediment to a wider use of SR systems (Lerman et al, 2004). Systematically varying individual agent parameters to study their effect on collective swarm behavior is impractical, even in simulation. Mathematical characterization of collective swarm behavior is more difficult, but provides the means to precisely predict it a priori-without the need of repetitive experiments.
Given the complexity of SR systems, and the frequently non-linear ways in which behaviors can arise, it is difficult to obtain precise bounds on the collective behavior of a swarm S of N robots each running a control algorithm κ. Robots might need to respond to environmental signals that arrive at unpredictable times; such systems are well-modeled as asynchronous, and therefore difficult to predict precisely. However, if we conceptualize S as a differentiable, continuous quantity, its dynamics can be modeled with Ordinary Differential Equations (ODEs) whose variables are the population counts associated with the different roles. We can apply a macroscopiccontinuous ODE modeling approach for the average behavior of S in the steady state (Berman et al, 2007), with the caveat that when using such a model to determine the behavior of S, actual system behaviors could be far from the average (Lerman et al, 2004). Usually, the larger the system, the smaller the fluctuations; in small systems the fluctuations resulting from finite size effects can be of order N, resulting in models which are accurate at asymptotically large scales but not at small scales. The master equation (Kampen, 2007), which is typically used to model the expected average behavior of systems, can be used to calculate the deviation from the average, but such calculations are often intractable or algebraically difficult.
A promising mathematically rigorous methodology utilizing macro-and microscopic models for group dynamics and individual behavior over time has been developed (Lerman and Galstyan, 2002;Berman et al, 2007;Galstyan et al, 2005;Sugawara and Sano, 1997), which sidesteps the difficulty in modeling the average behavior of the swarm, by instead modeling the change in the average behavior of the swarm, which is much easier. It uses (1) differential equations to model the behavior of the average number of robots in S in a given state, (2) discrete difference equations to model the stochastic transitions between robot states, and (3) stochastic simulation of discrete difference equations to compute state transition rates for all robots. It draws on implementations of the stochastic master equation in chemistry and statistical physics (Kampen, 2007). Through the usage of rate constants and population fractions in each state, it is possible to assess the general behavior of a system under a variety of stochastic circumstances. Most importantly, this approach has predictive control and performance guarantees-crucial components to translating laboratory models into viable real-world solutions without needing behavioral characterization of robot swarms via simulation experiments. Models in this differential equation paradigm operate on both the forward problem, i.e., predicting collective behavior from features of the control algorithm each robot runs (Lerman and Galstyan, 2002), and the inverse problem, i.e., incorporating design constraints into algorithm design in order to produce a desired collective behavior (Berman et al, 2007;Hsieh et al, 2008). However, the models are not overly robust, for two reasons. First, they make large simplifying assumptions such as homogeneous agent distributions, homogeneous environments (e.g., no obstacles and/or a completely visible arena), and Markov/semi-Markov scenarios (Berman et al, 2007). Second, they require many free parameters and extensive post-hoc model fitting which are specific to the implementation of the algorithm under study. Nevertheless, many notable applications of various forms of this methodology have appeared in the literature, demonstrating its practical utility. A few noteworthy examples include the stick pulling experiment (Ijspeert et al, 2001), foraging of green/red pucks using agent memory (Lerman and Galstyan, 2003), the "house hunting model" (Hsieh et al, 2008;Berman et al, 2007), and ant-inspired models that collaborate with or without communication (Sugawara and Sano, 1997).

Background in ODE modeling and diffusion
To model swarm collective behavior using ODEs, we consider the individual robot Finite State Machine (FSM) (Fig. 1). Our FSM is identical to previous work (Lerman et al, 2001;Lerman and Galstyan, 2002). Each state maps directly to a single robot behavior or a set/sequence of robot behaviors which together make up the robot controller for executing a foraging task. It is a coarse-grained model of robot behavior, which omits controller details such as sensing and actuation and contains the minimum number of states needed to describe the system dynamics for the problem. Such an approach is generally more mathematically tractable. If the results of the analysis do not sufficiently agree with the observed collective behavior, other states can be added (Lerman and Galstyan, 2002).  Another benefit of the ODE modeling approach is that each of the states in Fig. 1 maps directly to an ODE describing this change, where the transition rates for states become the ODE terms, and all ODEs can therefore be directly written down from robot control algorithm characteristics. We refer the reader to (Lerman et al, 2001;Kampen, 2007) for the theoretical underpinnings of the validity of this translation, which can be applied to any SR system, not just those doing a foraging task. The quantities modeled in Fig. 1 are listed in Table 1, and the ODE model of these quantities from previous work (Lerman et al, 2001;Lerman and Galstyan, 2002) is summarized in Eqs. (1) to (4). In our model, shown later in Section 4.3, we simplify the original equations and replace some parameters with mathematical derivations. Eq.
(1) describes the change in the number of robots in the searching state, which decreases as searching robots pickup blocks or encounter other robots and switch to collision avoidance.
The rate at which robots leave the searching state and switch to collision avoidance, α r N s (t)[N s (t) + N], can be understood as follows. When a searching robot encounters another searching robot, both switch to the collision avoidance state, decreasing the N s (t) by two; when a searching robot encounters either a homing robot or a robot already in the collision avoidance state, N s (t) decreases by one. The total decrease is then proportional to 2N s (t)+N h (t)+N s av (t)+N h av (t) = N s (t) + N, where N is the total number of robots in the swarm. Finally, searching robots encounter other robots at rate α r , regardless of the state of the other robot. Eq. (1) increases as homing robots deposit blocks in the nest or as searching robots exit the collision avoidance state. Eq.
(2) describes the change in the number of robots in the homing state, which increases as robots acquire and pickup blocks or leave the collision avoidance state, and decreases as robots enter the collision avoidance state or deposit their block in the nest.
can be understood analogously to the one in Eq.
(1). The model assumes that the rate at which homing robots encounter other robots is different than for searching robots, reasoning that there will be more congestion near the nest, and therefore uses a separate parameter α r to account for this. Eq.
(3) describes the change in the number of robots avoiding collision with other robots.
Eq. (4) describes the change in the number of free blocks available for robots to find; it decreases whenever a searching robot acquires a block (i.e., once deposited in the nest, blocks are not re-distributed in the arena).
In addition to using ODEs to model collective behavior, we also draw on diffusion theory. In diffusing systems composed of homogeneous particles undergoing normal diffusion, the average particle displacement is proportional to the diffusion time, i.e., a linear relationship; this assumes the particle is moving in an infinite, structureless medium close to equilibrium. However, there are important circumstances in which the relationship between average particle displacement and time is non-linear, referred to as anomalous diffusion (Oliveira et al, 2019;Metzler et al, 2014;Vlahos et al, 2008). For instance, in biological systems, interactions with other particles or membranes could influence the observed diffusion for macro-proteins as they traverse through biological media (Santamaria-Holek and Vainstein, 2009;Weiss et al, 2003;Nicolau Jr et al, 2007). Anomalous subdiffusion arises due to crowding in a concentrated system, which can make it heterogeneous and disordered (Surya K Ghosh and Metzler, 2016). Crowded environments can obstruct the path or make certain spaces inaccessible to the diffusing particle and prompt the rise of this behavior.

Changes to previous work
The ODE model described in Eqs.
(1) to (4) from (Lerman and Galstyan, 2002) has some limitations which we address in our model.

Uniformity assumptions
In previous work S is assumed to be uniformly distributed in 2D space, and blocks to be scattered randomly; therefore, the model is only accurate on scenarios meeting these criteria. In Fig. 2 we show some of those cases. Random (RN) block distributions (Sugawara and Sano, 1997; Hecker and Moses, 2015) ( Fig. 2d), are appropriate in scenarios such as order fulfillment in a warehouse, but many scenarios cannot be modeled by such distributions (Campo and Dorigo, 2007;Sugawara and Sano, 1997;Castello et al, 2016;Hecker and Moses, 2015;Lerman et al, 2001). For example, transferring material from one side of a building to another location requires a single source (SS) ( . For other scenarios, such as evacuation of civilians from a disaster zone, the block distribution cannot be inferred a priori, and a power law (PL) distribution ( Fig. 2c), in which blocks are clustered in groups of various sizes, is appropriate (Hecker and Moses, 2015;Harwell et al, 2020). By considering the distribution of blocks in a general way using spatial clusters, we guarantee that the model itself describes underlying characteristics of the environment and is robust enough to be applied to a wide range of foraging scenarios.

Invalid steady-state assumptions
The model assumes that S has reached steady-state, but at t = 0 assumes a finite number of objects for the swarm to find (Eq. (4)). This does not model accurately long-running SR systems or those with a small number of blocks relative to N. We formulate our ODE model to handle both steady-state and non steady-state behaviors by considering two types of operating environments. First, those in which any block deposited in the nest is re-distributed immediately in the operating area, so the overall number of blocks available for robots to find in the arena remains the same over time. Second, those in which any block deposited in the nest is not re-distributed, resulting in fewer blocks available for robots to find over time, which is the paradigm used in (Lerman and Galstyan, 2002). Our model treats non steady-state foraging environments as a special case of steady-state behaviors by setting to zero the number of areas where blocks can be distributed after nest deposition (see A j later in Eq. (8)). In this work, we are interested in steady-state behaviors.

Many free parameters
Free parameters do not invalidate the theoretical basis of the model, but limit its reuse by requiring iterative parameter refinement. The need to do posthoc model fitting requires experimental simulations, which makes it difficult to use the model as a predictive tool. We instead derive analytical expressions for all parameters except τ av from scenario characteristics, greatly reducing the need for post-hoc model fitting and reducing the number of free parameters from five to one. τ av is directly computable from the robot controller characteristics, or otherwise obtainable via experiments with a single robot. We introduce C ca (m), to characterize the inter-robot interference in scenario m. Although it is a free parameter, it has a strong semantic context, making fitted values more meaningful and suggesting the possibility of future derivation.

Preliminaries
To make accurate predictions of collective swarm behavior in the steady-state for each of the block distributions shown in Fig. 2, our model makes the following assumptions: 1. The swarm density ρ S is relatively low, so the behavior of N robots can be approximated by a linear function of the behavior of a single robot. An important consequence of this is that we can ignore the effect of a robot avoiding collision and encountering another robot during its avoidance maneuvers. This is not necessarily a limiting assumption, as many real-world applications require sparsely distributed swarms (Tarapore et al, 2020). We relax it as part of our model evaluation experiments in Section 6. 2. N is sufficiently large such that S can be approximated using the mean field model. We relax also this assumption as part of our model evaluation experiments Section 6. 3. S is homogeneous, and the control algorithm used by the robots is a Correlated Random Walk (CRW) (Renshaw and Henderson, 1981;Harwell et al, 2020), which is a random walk where the direction of the next step is biased based on the direction of the previous step. The bias angle θ is drawn from a probability distribution f (θ ) (Codling et al, 2010). Robots do a CRW at an average velocity v s until they acquire a block, which they then transport to the nest using phototaxis, i.e., motion in response to light, at a constant speed v h . Robots have no memory. 4. To reduce congestion in the nest, robots do not return to the exact nest center to drop their carried object, i.e., they employ a reactive avoidance strategy. Instead, they choose a random point along their trajectory from where they enter the nest to treat as the nest center. This shortens the homing distance, and can be computed as shown later in Section 4.4.
5. S has reached steady-state at some t 0. This maps naturally to swarms with long-running autonomy, such as those deployed on large-scale agriculture and industrial automation applications.
None of the assumptions we made reduces the utility of our model in its application to real-world problems; i.e., they are consistent with the constraints imposed by the problems themselves and commonly used robot hardware. We will relax some of them during the evaluation of our model, showing that while they are helpful during derivation, they are not essential for practical applications of our model.

ODE model
Let M = {SS,DS,PL,RN} be the set of scenarios based on the block distributions shown in Fig. 2. In each m ∈ M, let the area where blocks can be distributed be a subset of the overall area A of the arena. Let j = 1, . . . , J be the sub-areas within A in which blocks can be distributed, each described by a tuple (A j ,c j ,d j ,ρ j ). A j is the area occupied by the sub-area j, c j is the center of the sub-area, d j is its dimensions, and ρ j is the mean steady-state block density within the sub-area. The area within A where blocks can be found is the union of these disjoint subsets: The value of ρ j varies across sub-areas in our model, so that the block encounter rate α b can be captured accurately even in extreme non-homogeneous block distributions. As an example, consider a degenerate case in which there are two areas, both represented by a circular sector around the nest. Robots would likely forage only from the inner circle and rarely reach the outer circle, and so a single ρ j is not sufficient to capture collective dynamics in such cases. We require both ρ j and the robot spatial occupancy distribution in order to accurately compute α b in a general way (see Section 4.5).
(1) to (3) in previous work described earlier in Section 3.
Eq. (8) uses the described block modeling method to capture the underlying behavioral dynamics of the swarm. We made the following additional assumptions about the block distribution for its derivation. First, whenever a block is redistributed all j are selected with probability proportional to the fraction of distributable area they contain. Second, blocks are distributed uniformly within a given j. Third, every j can hold any number of blocks, allowing two blocks to occupy the same location (i.e., stacking). We note that under our-steady state assumption, Eq. (8) can be solved analytically given B(0) and does not need to be solved numerically. We also observe that from Eq. (8) the number of blocks in the arena available for robots to find can increase and decrease, addressing one of the weaknesses in the model from previous work.
Next we derive analytical models for τ h ,α b , and α r from the arena geometry, number of blocks, block distribution, etc. We do not derive τ av , because it depends intrinsically on the interference avoidance strategy employed by S and therefore cannot be derived independently from κ without additional assumptions.

Derivation of homing time τ h
We build the intuition behind the block acquisition probability for a robot at location x in the arena as follows, modeling the nest as a single point x n . Since searching begins from the nest, the density of N s (t) must be greater near the nest because robots perform biased random walks which originate from a common point. Consequently, this non-uniform swarm spatial distribution means that the mean distance from the nest at which a searching robot encounters a block is not the same as the mean distance of a block from the nest. From (Codling et al, 2010), we have that the spatial occupancy distribution from the central point, as a result of a biased random walk with bias distribution Uni f orm(−θ ,θ ), falls off linearly, and we would intuitively expect the following: 1. It is moderated by ρ j = B j (t)/A j , because the rate of decay of the mean block acquisition distance as a function of distance from the nest within a given A j is slower for low ρ j . Consequently, a j sparsely populated with blocks will have minimal effect on the overall swarm block acquisition probability distribution, while a j densely packed with blocks will create an area of higher acquisition probability. 2. ρ j would play an exponentially moderating role only when the block acquisition location x is close to x n , such as for RN or PL block distributions. For SS and DS distributions, where the mean distance from a block to the nest is large, the effect of ρ j on block acquisition locations should be minimal.
Using the results of (Codling et al, 2010) and these intuitions, we formulate Eq. (9) as a close approximation of the occupancy distribution of a single random walker performing a correlated random walk starting from a central point x n . C is a normalization constant to ensure p acq j (x) integrates to 1 over all j.
Having defined the probability density function, we now derive the expected acquisition location by finding the expected values of the marginal density functions in x by integrating Eq. (9) and summing over all j: and similarly for y. We now write an expression for τ 1 h : v h is the robot phototaxis velocity, specified in the input configuration, and E[x acq ] − x n is the expected distance an acquired block will be from the center of the nest. d cr is the distance the homing path is shortened due to the employed congestion reduction strategy, and is straightforward to calculate from the arena geometry. We provide an example of computing d cr for a square arena, which in this work corresponds to the RN and PL scenarios; similar calculations can be done for the other arena shapes. Using a square nest of length L, without loss of generality we find the average distance to the origin of a randomly selected point in a region R := {(x,y) : 0 ≤ y ≤ x ≤ L/2} that is a triangle 1/8 of the nest. The uniform density on this region is then: Then, the average distance d cr is computed via: Finally, to derive τ h , we note that under our assumption of low to moderate ρ S , the homing time increases linearly with N according to the expected value of time lost due to inter-robot interference (Lerman and Galstyan, 2002), averaged across all robots: 4.5 Derivation of block acquisition rate α b Using our mean-field assumption, S can be approximated as a fluid composed of robot particles, and be considered to obey many of the same laws; in the long-term limit, the governing equation for the biased random walk used in this work is the advectiondiffusion equation (Codling et al, 2010). Given sufficiently simple robots, this approximation gives good results (Codling et al, 2010;Pang et al, 2019). Using simple robots does not detract from the utility of the system, which mimics the self-organizing behavior of real swarm populations such as social insects. Using this intuition, we obtain α b by computing the mean time it takes a robot "particle" starting at x n to "diffuse" within the operating area A to the expected acquisition location E[x acq ]. Viewing E[x acq ] − x n as the Root Mean Square (RMS) displacement distance and assuming a linear relationship between displacement and diffusion time, we obtain: where F(N) is the diffusion constant for a swarm of N robots, and α b is the expected time to diffuse from x n to E[x acq ]. We note that while the calculation for α b does not directly consider the distribution of blocks in the arena and their densities within a given block cluster, it depends on E[x acq ], which does.
An exact calculation of F(N) is out of the scope of this paper, so we approximate it as shown in Eq. (16) using the results of the RMS diffusion for CRW (Codling et al, 2010) and our intuition that smaller θ will result in quicker swarm diffusion (more straight line motion): where C d f (m) is a per-scenario parameter characterizing the linearity of the diffusion rate, and D xy is defined as follows from (Codling et al, 2010): where v s is the robot searching velocity and a scenario parameter. Under "normal" circumstances with inert particles, Eqs. (15) and (16) give good results because system diffusion varies linearly with time, which these equations require. However, many foraging environments are heterogeneous and have non-uniform distributions of blocks. As a result, the swarm may experience non-linear anomalous diffusion (Hasnain et al, 2018;Woringer et al, 2020), i.e., crowding in some areas of the environment, giving rise to the need for C d f (m), which implicitly captures these artifacts, as shown in Eq. (16).
To calculate C d f (m) from first principles, we take the following approach. First, we calculate the environment heterogeneity H(m) in reference to an environment with perfectly uniform block distribution (homogeneous case). For the square RN and PL scenarios, we compute the variance of the average inter-cluster distance. For the rectangular SS and DS scenarios with J ≤ 2, the average inter-cluster distance is not meaningful, so we calculate H(m) as the mean cluster distance from the "ideal", which is defined as an environment in which all c j are incident with x n . Modulating H(m) in Eqs. (18) and (19) by the expected level of crowdedness within each cluster we obtain an expression for C d f (m) which accounts for the number of block clusters (J) and the area cumulatively covered by all clusters (A d ): In Fig. 3, we plot the calculated values for C d f (m) for all scenarios to verify we have restored linearity of diffusion, and see strong correlations between C d f (m) and F(N) for all block distributions. As expected we see the highest heterogeneity associated with the lowest diffusion rate (PL), and the lowest heterogeneity associated with the highest diffusion rate (RN). With more observed heterogeneity within the environment, the scaling factor necessary to restore the linearity of the diffusion rate increases to accommodate for this occurrence. We also note that the differing x intercepts depend on the block distribution, and that as the size of the arena increases the crowding observed is not as pronounced (less anomalous diffusion). This might be due to the decreased likelihood of robots being "trapped" (Chepizhko and Peruani, 2013) by obstacles in the environment, leaving the robot more able to continue moving through the arena.

Derivation of robot encounter rate α r
To derive α r , we use our assumption of low to moderate ρ S , to model the robot encounter rate as a function of α 1 r and F(N). We view Fig. 1 as a queueing network, where robots are either performing collision avoidance maneuvers or not. The input/output transition rates for a state are summed to form the arrival and service rates for the collision avoidance queue Q ca (λ ca = α 1 r and µ ca = 1/τ av , respectively). Modeling Q ca as a M/M/1 queue, i.e., at most one robot exits collision avoidance per ∆t, which is reasonable if ∆t is small, we can write the following using Little's Law (Šeda et al, 2017) as: The second term in Eq. (20) is a corrective factor accounting for robots that experience interference due to encountering arena walls, not other robots, which is simply the scaled rate at which a single robot experiences interference. We cannot use N av (t) directly, as α r needs to be computed a priori, but we can estimate it as a function of α 1 r and N 1 av (t), using our intuition regarding swarm diffusion: where α 1 r and N 1 av (t) can be computed from κ using the results of (Codling et al, 2010). We increase the influence of θ in Eq. (21) by introducing another D θ in the denominator; smaller θ will result in more interrobot interference due to straight line motion. C ca (m) characterizes the sub-or super-linearity of Eq. (21) for a scenario m in relation to the random (RN) scenario.

Experimental setup
We use the ARGoS simulator (Pinciroli et al, 2012) 1 with a dynamical physics model of the marXbot robot in a 3D space for maximum fidelity (robots are still restricted to motion in the XY plane).
Robots run the FSM shown in Fig. 1, use short range proximity sensors for collision avoidance, and ground sensors to check if they are currently on top of a block. Blocks are abstracted as black squares on the arena floor. We ran four sets of experiments characterized by different swarm density, and evaluate our model across the scenarios shown in Fig. 2, relaxing some core assumptions made when deriving our model as we proceed.
We ran our experiments with N = 1 . . . 6,006 robots in environments with SS and DS source block distributions, and with N = 1 . . . 11,837 robots for RN and PL block distributions. Our scenarios are obstaclefree in the sense that they do not contain external obstacles; however, robots act as obstacles to each other. We use θ = π 36 for CRW and compute τ av from the details of our robot control algorithm κ (details omitted). C ca (m) is computed post-hoc for each experiment. For all experiments we performed 32 runs of T = 200,000 seconds, with a timestep of t = 0.2 seconds, i.e., the robot control algorithm runs at 5 Hz.

Results
For all the scenarios, we omit predictions for N s (t) since it can be computed from the conservation of robots via N s (t) = N − N av (t) − N h (t). We present our results with a non-scaled X-axis, with one data point per experiment, and a logarithmic Y-axis for larger N, to improve readability.
6.1 Constant ρ S , large N These are the "ideal" conditions under which to evaluate our model, because the linearity and mean-field assumptions we made when deriving our model hold.
In Fig. 4, we see strong agreement between the predictions of our model and experimental results for all scenarios, providing compelling evidence that our model is capturing the underlying dynamics of interrobot interference and searching accurately and that our mean-field and linearity assumptions are accurate. PL scenarios are the least favorable of all foraging environments, since they are asymmetrical and do not contain easily exploitable block clusters. Our model struggles to predict N h (t), N av (t) within the 95% confidence interval for PL scenarios, but does track the overall trend reasonably well, showing that our underlying diffusion model and assumptions about linearity of α r are generally accurate even in this difficult case. The divergence that does exist between predictions and experiments for PL scenarios suggests that Eq. (9) is moderately inaccurate; this is further supported by slight differences between experimental data and predictions for N h (t) in the RN and SS scenarios (a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for RN and PL scenarios. environments. In Fig. 4b, when our predictions of N av (t) are inaccurate for some DS scenarios (in terms of confidence intervals), the inaccuracy may not be of practical concern as the predictions only differ from experimental results by < 1 robot/6,000.

Constant ρ S , small N
We relax the mean-field assumption, and evaluate our model's ability to capture finite-size effects when S cannot be reasonably approximated using the meanfield model (with N ≤ 50).
In Fig. 5, we see strong agreement between experimental results and the predictions of our model for all tested scenarios with small N. Even when this crucial assumption is relaxed, our model accurately captures finite-size effects, suggesting that the mean-field assumption may be valid at very small N, depending on the control algorithm used. The greatest discrepancy between predictions and results is shown in Fig. 5c for RN scenarios, which are inherently the most non-linear because of the placement of the nest relative to the objects to be gathered. In comparison (a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) with experimental data for RN and PL scenarios. with results from the previous section, stronger agreement between predictions and results with constant ρ S is seen when N is small, likely due to the inherent stochasticity of swarms which comes into play more strongly at larger N.

Variable ρ S , small N
We relax the linearity and mean-field assumptions, and evaluate our model's ability to capture finite-size effects when S exhibits non-linear behaviors with ρ S = 0.01 − 0.1 and N ≤ 50.
In Fig. 6 we see that our model generally performs well even when we violate its low ρ S assumption. We see strong agreement between the model and experimental data for all scenarios for ρ S = 0.01 − 0.05 for both N h (t) and N av (t). In Fig. 6d we see a striking divergence between our model and results at ρ S = 0.09, which may be a numerical anomaly or the point at which our ODE solver begins to struggle with the non-linear behavioral inputs.
(a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) N av (t) = N h av (t) + N s av (t) (avoidance) with experimental data for RN and PL scenarios.

Variable ρ S , large N
We restore the mean-field assumption to evaluate our model's ability to capture collective dynamics across scales on the order of 1,000s of robots when S exhibits non-linear behaviors. We believe that ρ S = 0.1 is an appropriate upper limit for real-world swarms of practical utility, as it corresponds to 1 robot/10m 2 which is an extremely high density in the context of feasible real-world applications.
In Fig. 7, we see that our model generally does not model swarm behavior at large scales well. In Fig. 7 we see divergence between model predictions and results for N av (t), N h (t) for all scenarios beyond ρ S = 0.02, and beyond ρ S = 0.01 for RN, PL scenarios, and our model fails to provide meaningful predictions, due to the highly non-linear swarm behaviors present, which is as we expected. For RN and PL scenarios, we see a strong correlation between the ρ S at which our model's predictions for both N h (t) and N av (t) break down, which happens at ≈ ρ S = 0.01.
(a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for RN and PL scenarios.

Discussion
Our results advance the state of the art in ODE modeling for foraging swarms (Lerman et al, 2001;Lerman and Galstyan, 2003) by demonstrating accurate predictions of swarm behaviors across scales and block distributions when the assumptions made in Section 4.2 are valid. To the best of our knowledge, our model is the first presented in the literature capable of capturing the collective dynamics of swarms when the mean-field approximation is valid, i.e., large N, but also when the mean-field approximation is not valid and finite-size effects are present. i.e., small N. Building on our results, accurate predictions of swarm performance and SR system properties, such as scalability and emergent self-organization using (Harwell et al, 2020), are possible.
As an example, we can predict swarm performance with reasonable accuracy using our model, as shown in Fig. 8. We define a performance measure P as the (a) SS and DS: N = 1 . . . 50.
(b) RN and PL: N = 1 . . . 50. Fig. 8 Predictions of swarm performance from first principles for different scenarios rate of block collection by the swarm S. We then predict it via P = α b /µ h , where µ h is the rate of robots leaving the homing queue Q h (analogous to Q ca ). Under our assumption of Q h as a M/M/1 queue, we set µ h = 1. In Fig. 8 we see strong agreement between our predictions and experimental results for SS and PL scenarios across all swarm sizes, and reasonable agreement for DS scenarios. Our predictions are less accurate for RN scenarios, likely due to our approximation of Q ca as an M/M/1 queue not being correct; i.e., our model does not account for multiple robots dropping an object in the nest on a single t. When we consider the swarm behavior with constant ρ S at large N (Fig. 4), we see our model struggle to predict behavior for all portions of the curve for all scenarios. We see that while our model predicts behavior at N ≥ 100 well, it struggles to predict behavior at N < 100 for SS, DS, and RN scenarios, underscoring the difficulty in creating a single linear model which can capture both mean-field and finite-size behaviors simultaneously. PL scenarios are the most challenging to model, and we see that while our model follows the general trend of behavior well, its predictions are usually outside of the confidence interval range. However, our model is accurate to within ≈ 16 robots/11,837, which may be sufficient in many cases.
Our linearity assumption, strictly speaking, is clearly not valid in all contexts across scales or across scenarios at similar N. Nevertheless, it is valid in many contexts (e.g., Figs. 5 and 6, parts of Fig. 4), and thus our model is of practical utility in designing SR systems to tackle problems across scales. Furthermore, our findings suggest that N always has a non-linear effect on behavior, even at extremely low densities, which is a function of the scenario itself. For more favorable scenarios, this non-linearity may be negligible even at large N, but for less favorable scenarios, such as RN and PL, it is not, and non-linear modeling must be applied for accurate prediction across scales.
When we relax our linearity assumption and increase ρ S to introduce non-linear behaviors in the swarm via inter-robot interference, our model still produces comprehensively accurate predictions of behavior at small N for all scenarios (Fig. 6). While our model will undoubtedly break down at some ρ S > 0.1, our results demonstrate that the threshold is above any realistic scenario, because densities > 0.1 are infeasible in general for real-world problems.
Capturing non-linear behaviors using linear models is much more successful if N is small; this aligns well with our intuition. If we imagine that the amount of time each robot spends avoiding collision with other robots at each timestep (t ca ) is roughly of the form t ca = N ρ S then it is clear that t ca will increase more quickly with N, even if ρ S is low.
Based on this discussion, we can confidently characterize the limits of linear modeling, at least for ODE based approaches, and say that it is possible to capture non-linear behaviors of swarms using linear models across scales. Furthermore, for N ≤ 50, doing so is also practical, as our model still provided predictions accurate within 95% confidence in all scenarios. However, scenarios exhibiting comparatively more nonlinear behaviors via N av (t), such as RN, have large 95% confidence intervals, reflecting the greater variability in behavior as a result of non-linearity. As a result, our model's predictions of the swarm average behavior may be less useful in such scenarios, since the model provides predictions of the average case. We note that while the results of this work suggest that linear modeling of non-linear swarm behaviors is possible using other modeling approaches, they do not necessarily generalize to other modeling paradigms.
Our model still relies on one parameter: C ca (m). While this is a free parameter in the strictest sense, it has a strong semantic context rather than just being a fitting parameter to ODE terms to improve prediction accuracy. This is a substantial improvement over previous work which used five context-free fitting parameters (Lerman and Galstyan, 2002). Thus, while our model is not (yet) strictly fundamental, the semantic context of C ca (m) strongly suggests that it can be determined analytically, and first-principle modeling of general foraging behavior is possible.
Finally, it is important to emphasize the presented results were obtained in simulation, commensurate with other mathematical approaches to modeling SR systems (Guerrero-Bonilla et al, 2020;Tarapore et al, 2017). Even if our model does not produce equivalently accurate results to those shown in Figs. 5 and 6 with systems of real robots, we believe that it will still provide a strong starting point for modeling systems of real robots, thereby accelerating the deployment of swarm solutions to real-world problems, for two reasons. First, the robot FSM used in this work (Fig. 1) is simple: it has minimal sensor usage, and does not use memory. This simplicity will make real-world systems of such robots especially amenable to mathematical modeling. Second, our model has been shown to be accurate with swarms of non-trivial size (e.g., up to 46 robots in Fig. 5, up to 51 robots in Fig. 6, parts of Fig. 4 with up to 11,837 robots). In such swarms, nondeterministic transient behaviors will reliably arise even in simulation, e.g., floating point representation errors; similar errors will arise even in small swarms of real robots.

Conclusions and future work
We have investigated how to characterize the limits of linear modeling of swarm behaviors. We have developed a robust ODE model for foraging swarms, and shown it is accurate in a wide variety of scenarios in swarms which exhibit linear behaviors, and are large enough to be approximated using mean-fields, demonstrating its utility as an effective modeling tool. Furthermore, in swarms for which either assumption is not true, our model in many cases is still accurate. Future work will derive the block acquisition density function using random walks theory and refine our model with experiments using real robots.