Decentralized self-selection of swarm trajectories: from dynamical systems theory to robotic implementation

In this paper, we present a distributed control strategy, enabling agents to converge onto and travel along a consensually selected curve among a class of closed planar curves. Individual agents identify the number of neighbors within a finite circular sensing range and obtain information from their neighbors through local communication. The information is then processed to update the control parameters and force the swarm to converge onto and circulate along the aforementioned planar curve. The proposed mathematical framework is based on stochastic differential equations driven by white Gaussian noise (diffusion processes). Using this framework, there is maximum probability that the swarm dynamics will be driven toward the consensual closed planar curve. In the simplest configuration where a circular consensual curve is obtained, we are able to derive an analytical expression that relates the radius of the circular formation to the agent’s interaction range. Such an intimate relation is also illustrated numerically for more general curves. The agent-based control strategy is then translated into a distributed Braitenberg-inspired one. The proposed robotic control strategy is then validated by numerical simulations and by implementation on an actual robotic swarm. It can be used in applications that involve large numbers of locally interacting agents, such as traffic control, deployment of communication networks in hostile environments, or environmental monitoring.


Introduction
How can one trigger the self-emergence of desired spatio-temporal patterns via a collection of interacting autonomous agents?Can we analytically unveil such collective behavior for scalable collections of agents?Is it possible to explicitly demonstrate how such collective behavior can possibly remain stable in the face of environmental noise?Progress in the engineering of swarm robotics relies on explicit and detailed answers to these questions, and that is the central goal of the present contribution.As no generic guidelines yet exist, engineering such collective behavior remains, in many respects, an open challenge.Bonabeau, Dorigo and Theraulaz illustrate how elementary, fully decentralized biological systems can achieve highly elaborate collective tasks Bonabeau et al. (1999).This pioneering work remains a rich and challenging source of inspiration for many further developments both in ethology [see the recent overview by Leonard (2013)] and in swarm robotics.Presently, serious effort is being made to import some basic features of these complex living systems into the engineering world of robotics.Obviously, animals remain intrinsically highly complex machines compared to actual robots and, therefore, direct applications remain rather limited.Nevertheless, recent attempts to combine artificial and natural collective systems by implementing the models observed in animal societies into that of robots shows promise (Halloy et al. 2007;Mondada et al. 2011).Because of the complexity of the interactions taking place in such mixed societies, the models often remain partial and the link between them has not been fully established.
Several up-to-date reviews (Kernbach 2013;Brambilla et al. 2013) summarize ongoing research in this area and the various engineering applications.Often the design of such robotic systems consists of a bottom-up trial-and-error exercise with no analytical link between the microscopic model (describing the single agent behavior) and the macroscopic one (describing the global result achieved).Microscopic and macroscopic models are then seen as separate approaches Brambilla et al. (2013), rather than steps in the engineering process.To bridge the gap between the micro-and macroscopic descriptions, Martinoli et al. (1999) constructed the first link between individual behavior and the statistical macroscopic model, which was successfully implemented on real robots.Along the same lines, Lerman et al. (2005) exhibited the emergence of collective behavior from individual robot controllers using a class of stochastic Markovian mathematical models.The authors validated their approach by performing experiments using real robots.More recently, Hamann and Wörn (2008) used an explicit representation of space together with methods of statistical physics to establish the link between microscopic and macroscopic models.Space heterogeneities are also considered by Prorok et al. (2011), who derived, from diffusion-type evolutions, a collective inspection scenario implemented on real robots.Brambilla et al. (2012) proposed a top-down design process built on iterative model construction and based on Probabilistic Computation Tree Logic (Deterministic Time Markov Chains).Their design methods were further validated on a group of physical e-puck robots; their iterative construction potentially enables other types of swarming behavior, depending on the skills of the acting programmer.Berman et al. (2009) likewise proposed a top-down design methodology, but at a higher level (i.e., at the sub-swarm level).Their decentralized strategy allows dynamical allocation of a swarm of robots to multiple predefined tasks.The approach by Berman et al. (2009) is based on population fractions, which enable the design of stochastic control policies for each robot, and does not assume communication among robots.This strategy is analytically proven to achieve the desired allocation, and validation is made on a surveillance scenario involving 250 simulated robots distributed across four sites.More recently, Berman et al. (2011) have extended their work to spatially heterogeneous swarms.This control, however, can only generate static patterns of robot sub-swarms.In the present work, we aim to generate self-organized dynamic patterns (i.e, circulation of robots along closed planar curves).
Closely related to the present paper, Hsieh et al. (2007) synthesized controllers that allowed individual robots, when assembled in a swarm, to self-organize and circulate along a predefined closed curve.The system is fully distributed and relies only on local information, thus ensuring scalability.The controller can be analytically described, and does not require communication between the robots, thus simplifying its implementation.Hsieh et al. (2007) demonstrated convergence of the system for a certain class of curves and validated their theory with simulations.While the final behavior looks similar to what we are targeting here, we aim at directly deriving our controller from a stochastic dynamical model that can be analytically discussed.Our goal here is to allow the swarm to converge to a consensually selected curve among a given family (and not to a predefined curve as done in Hsieh et al.'s study).Conceptually, Hsieh et al.'s study belongs to the important family of gradient-based controls.Within this family, Pimenta et al. (2013) proposed a decentralized hydrodynamics-inspired controller that drives the swarm into a preassigned arrangement using a static and/or dynamic obstacle avoidance mechanism.They also tested their approach on swarms of actual robots, including e-puck robots.More closely related to our paper, Sabattini et al. (2012) developed a control strategy using decentralized time-dependent potential construction.Specifically, this mechanism allows a robotics swarm to track a closed curve (given by an implicit function, as in this paper), while keeping a minimal safety distance between the robots.In Sabattini et al.'s study Sabattini et al. (2012), the robots are given a predefined curve, whereas in our approach, the swarm consensually selects the closed curve that will be tracked.
Among non-gradient-based approaches, Schwager et al. (2011) demonstrated the stability and quantified the convergence rate of the control of a swarm of quad-rotors.Their approach is based on a combination of decentralized controllers and mesh networking protocols.Moreover, the authors only used local information and validated their control on a group of three real robots.Schwager et al.'s controller needs to know the relative position of its neighbors with high precision, while our goal is to be noise-resistant.Sepulchre et al. (2007) developed a methodology for the generation of stable circular motion patterns based on stabilizing feedbacks derived from Lyapunov functions.This theoretical work was then used to design the control of mobile sensor networks for optimal data collection in environmental monitoring Leonard et al. (2007).The global patterns generated in Sepulchre et al's study exhibit similarities with those presented in the present paper, but their control is "less" decentralized.Namely, Sepulchre et al. assume an all-to-all communication network between the robots, whereas in our work, we aim at using only local communication between neighbors within a given range.Given Sepulchre et al.'s assumption, the controller is also more sensitive to losses of communication between robots in the swarm.Finally, within non-gradient-based controls, it is worthwhile to mention Tsiotras and Castro (2011), who devised a decentralized mechanism, leading to Spirograph-like planar swarm trajectories.Their approach only relies on individually gathered information on distance and angle to neighbors.
With respect to the state of the art, whereas most research in control theory aims at finding strong theoretical guarantees, the present work follows a very different approach.We aim at bridging the gap between the purely mathematical considerations and the implementation on actual robots, and approach the problem from a mathematical angle.The mathematical models are then implemented using a swarm of mobile robots.To some extent, our work is closer to that of biologists who use robotics as a validation tool for their theoretical models.Using this approach, we are able to achieve a decentralized consensual swarm control mechanism using only simple local information, such as the number of neighbors in a given range and their relative positions.In order to better link the model and the implementation, this article also includes a dictionary to help translation of the main concepts between the theoretical and the robotic control mechanisms.The developed controllers are not designed to achieve a fixed shape or trajectory, as often found in the state of the art, but to collectively choose one consensual shape.By tuning the control parameters, one can enable a given swarm to converge on different orbits without changing the behavior of each agent.Another important feature is the fact that the consensual orbit is not known by the agents, but is only determined by mutual interactions in a fixed range.
Although perfectly predictable analytically, this collective choice can be described as a collective weakly emergent behavior.Weak emergence is here understood in the S. and O. Kernbach's interpretation (Kernbach 2013;Kernbach et al. 2013) (i.e., macroscopic swarm behavior due to interactions between agents but having its origin in the model that describes the system's behavior).Our model has a tractable complexity, enabling us to explicitly describe the resulting emergent collective behavior, which is controllable.Specifically, our class of models presents a distributed mechanism, enabling agents to select a parameter while fixing one consensual traveling orbit among the given family in a distributed manner, and this feature holds even if noise corrupts their movements.
The paper is organized as follows: Sect. 2 exposes the basic mathematical formulation of agents' dynamics.Section 3 reports simulation experiments that fully match our theoretical finding even for finite N agents' populations (i.e., when the mean-field limit is not reached).Section 4 describes, in detail, how the actual implementation of our mathematical model was realized with a collection of e-puck autonomous robots.

Mathematical modeling
Our mathematical model consists of a set of N coupled stochastic differential equations (SDE) where the driving stochastic processes are White Gaussian Noise (WGN) representing a swarm of homogeneous Brownian agents Schweitzer (2003).The use of WGN in the dynamics implies that the resulting stochastic processes are Markovian.Accordingly, all probabilistic features of the model are known by solving the associated Fokker-Planck (FP) equation.A similar approach has already been used in several studies (Berman et al. 2011;Hamann and Wörn 2008;Prorok et al. 2011).Agents' interactions are ruled by specifying the radius of an observation disk (i.e., metric interactions) and the strength of the attraction/repulsion control parameter.The specific form of agents' individual dynamics and their mutual interactions are constructed from a Hamiltonian function, from which one derives both a canonical and a dissipative drift force (i.e., mixed canonical-dissipative vector fields Hongler and Ryter (1978)).By imposing the orthogonality between the canonical and dissipative vector fields, we are able, in the N → ∞ mean-field limit, to explicitly calculate the stationary probability (i.e., invariant probability measure) characterizing the global agents' populations.For our class of models, we can analytically observe how a decentralized algorithm is able to generate a milling-type spatio-temporal pattern.In this pattern, all agents will circulate with constant angular velocity in the vicinity of a self-selected level curve from the input Hamiltonian function.

Single agent dynamics
Let us first consider the single agent Mixed Canonical-Dissipative (MCD) stochastic planar dynamics (a tutorial devoted to MCD can be found in Appendix section "MCD dynamics"): where V stands for an antisymmetric (2 × 2) dynamical matrix defined by: where When σ = 0, the (deterministic) dynamics Eq. ( 1) admits closed trajectories C 's, which are defined by the solutions of the equation Depending on the sign of the curvature (i.e., second derivative ) U taken at C , one can either have attracting (U < 0) or repulsive (U > 0) limit orbits.An attracting C drives any initial conditions (X 0 , Y 0 ) to C in asymptotic time.In Eq. ( 1), the antisymmetric part of V produces a rotation in R 2 with angular velocity ω.This defines a canonical (conservative) part of the vector field (i.e., the Hamiltonian part of the dynamics).Concurrently, the diagonal elements of V model the non-conservative components of the dynamics (i.e., dissipation or supply of energy).From its structure, the non-conservative (deterministic) vector field of Eq. ( 1) is systematically perpendicular to the conservative one.This geometrical constraint reduces the dimensionality of the dynamics and allows the derivation of explicit and fully analytical results.
In the presence of WGN (σ > 0 in Eq. ( 1)), the resulting SDE describe a Markovian diffusion process, fully characterized by a bivariate transition probability density (TPD) measure where δ(•) stands here for the Dirac probability mass and the FP operator reads: For asymptotic time, the diffusion processes reach a stationary regime: lim t→∞ P(x, y, t | x 0 , y 0 , 0) = P s (x, y), solving the stationary FP Eq. ( 3), for which we obtain: and N is the probability normalization factor yielding R 2 P s (x, y)dxdy = 1.
Example Consider the Hamiltonian function In this case, the deterministic dynamics coincides with the Hopf oscillator with a unique attracting C that is a circle with radius L .In this situation, the measure Eq. ( 5) exhibits the shape of a circular "ring" with its maximum located on the circle C (see Fig. 1).

Multi-agent dynamics
Let us now build on the dynamics presented in Sect.2.1 and consider a homogeneous swarm of N interacting agents exhibiting individual dynamics similar to the one described in Eq.
(1).Specifically, for k = 1, 2, .., N , we now write: where V k,ρ,M stands for an antisymmetric (2 × 2) dynamical matrix defined by: where M > 1 and ω > 0 are (time-independent) control parameters, M •ω is a normalization factor and dW (•) (t) are all independent standard WGN sources (homogeneity follows as we ) is located at the center of an (individual) observation disk D k,ρ (t) with common constant radius ρ.We further assume that A k is, in real time, able to count the number of neighboring fellows N k,ρ (t) that are contained in D k,ρ (t).Finally, we introduce the instantaneous (geometric) inertial moment of the agents located in D k,ρ (t), namely: To heuristically understand the collective behavior emerging from the set of SDE's given in Eq. ( 6), let us comment on the following features: (a) Large population of indiscernible agents.As we consider homogeneous populations, we may randomly tag one agent, say T , and consider the actions of the others as an effective external field (often referred to as the mean-field point approach).This procedure reduces the nominal multi-agent dynamics to a single effective agent system, thus making the analytical discussion much easier while still taking into account interactions of the neighbors.The individual T -diffusive dynamics on R 2 are then assumed to be representative of the global swarm and follow from dropping the index k in Eqs. ( 6) and ( 7).The effective interactions terms N ρ,M and L ρ,M are to be determined by auto-consistency from the TPD measure of T .In the stationary regime, thanks to the orthogonal property between the canonical and the dissipative components of the vector fields, we obtain: , the harmonic oscillator Hamiltonian, and for large signal-tonoise ratio S = 2γ /σ 2 1, we can analytically express the value of the limit cycle's orbit: (b) Roles of the control parameters M and ρ.The origin O := (0, 0) is a singular point of the deterministic dynamics and its stability follows from linearizing the dynamics near O.In the stationary regime, the linearization reads: with eigenvalues: When RE < 0, (resp.> 0), O is attractive (resp.repulsive).The very existence of a stationary regime implies that P s (x, y) vanishes for large values of |x| and |y|.For agents far from O, both N k,ρ,M and L k,ρ,M will vanish (with high probability, except the agent at the center, no agents will be found inside D k,ρ ), implying O to be attractive as RE < 0; conversely, agents close to O will be repelled.Ultimately, an equilibrium The black circle (hidden below the swarm of agents) gives the theoretically computed limit cycle.This shows the exactitude of Eq. ( 10), even for finite swarms.Left Harmonic oscillator Hamiltonian function Yates (1947)).Right Same Hamiltonian function as the center figure, with ρ = 0.8, presenting one connected orbit around both attractors consensus will be reached when the diagonal part of V ρ,M vanishes which is achieved when Eq. ( 10) is attained.(c) Role of the dissipative terms U [H, t] and the parameter L ρ,M .Observe that the regulating term U [H, t] is itself not necessary to reach a consensual behavior.However, for large γ 's, it enhances the convergence rate toward equilibrium and reduces the variance around the mode of P s (x, y).This enables us to easily estimate the radius L ρ,M given by Eq. ( 10).
According to Eqs. ( 9) and ( 10), we observe that, in the stationary regime, the swarm selfselects a consensual circular orbit with radius L ρ,M that depends on both the observation range ρ and the effective attracting range M.
Numerical simulation results, obtained by applying the Euler method over discrete timesteps with dt = 5 • 10 −3 [s], are shown in Figs. 2 and 3.The employed discretization works perfectly, as the numerical results coincide very closely with our theoretical ones (especially in the case of the theoretical limit cycle in Fig. 2, left).
Note that a Hamiltonian function can be constructed from a set of points along any desired closed curve on the plane (without self-intersection).A simple fit can be made on this set of points (e.g., least squares method), returning a high-degree polynomial form representing the Hamiltonian function.This function will be differentiable (by construction), and we can ensure H (X, Y ) > 0(∀X, Y ) by adding a constant; the Hamiltonian can then be directly used in our model.

Braitenberg control mechanism
The actual adaptation of our MCD dynamics is implemented on the model of a Braitenberg control mechanism (BCM) Braitenberg (1984).In the BCM, the speed of the robot's motor(s) is reactuated at discrete timesteps, solely based on the output value(s) of a series of sensors.
In our case, we use cylindrical robots, equipped with eight light sensors (s0,s1,...,s7) evenly distributed along their circular perimeter and two driving motors (left and right), as pictured in Fig. 4. The robot's movement evolves on the plane, with a light source located at the origin.
, resulting in two disconnected orbits The light source serves as the origin O in the mathematical model, acting as an attractor for the robots, while also being the center of the limit cycle the robots will ultimately arrange and circulate on.Our goal is to let a swarm of robots arrange and circulate in a circle around the light source, following a similar control mechanism to Eq. ( 1), with the Hamiltonian function of the Hopf oscillator.This choice of circular trajectories allows us to compare the behavior of an agent swarm (with analytical results) and a robotic swarm.
The gathered information enables us to update the velocities at discrete timesteps according to the rule: where α and β are two positive parameters, controlling the translation/rotation velocities.This pair of controls in Eq. ( 12) effectively adapts Eq. ( 1) (in Cartesian coordinates) to a robotics paradigm, in terms of left/right motor speeds (speed and heading).In Eq. ( 12), we reduce the number of sensors to two by grouping them into the left group (S L ) and the right group (S R ).
The first term α • N ρ (t) accounts for the forward speed of the robot, and the other term regulates the heading.Indeed, when α = 0, the speeds have the same norm but opposite signs, implying that the robot can only turn around its main axis.In this case, the robot stabilizes with the light source facing S, i.e., with the light source at the same distance from sensors 0 and 1, with S L (t) ≡ S R (t) ⇒ V L (t) = V R (t) = 0.This ending position is the only stable steady state; when the light source is facing the point U , at the same distance from sensors 4 and 5, it is unstable.
When the light initially stands on the left of the robot (i.e., facing sensors 5, 6, 7, or 0), V L (t) will exceed V R (t) (as S L (t) > S R (t)).In this configuration, the robot will turn on its axis until the light faces the point S. The same argument holds when the light initially stands on the right of the robot.
When the first term in Eq. ( 12) (the forward speed) is constant and positive, this control mechanism drives the robot to the light source to turn around it.The robot constantly moves forward, keeping the light source facing S, thereby creating a rotating trajectory around the light source.
The radius of the resulting orbit will be directly proportional to α, and inversely proportional to β: an increase in α results in a larger distance covered before the next heading update (and hence, a larger orbit).Conversely, an increase in β induces a sharper heading adaptation (and hence, a smaller orbit).For a fixed β, the radius of the limit cycle depends only on α.
In our case, we fix α and β, thus letting the forward speed be proportional to the number of instantaneous neighbors.This forces the orbit to shrink when few neighbors are detected, and increase in the other configuration, mimicking the MCD dynamics.While in the mathematical model of Sect.2, the parameters ρ and M selected the limit cycle's radius, here ρ, α, and β contribute to the orbit selection.Accordingly, the parameter M will be adjusted by an ad-hoc combination of α and β.
Assuming the swarm of robots is large enough to ensure local communications, the proposed BCM will let the swarm converge to a consensual orbit and circulate on it at a constant speed.This convergence does not depend on the initial position and heading of the robots.As the central light is a global attractor, robots will also never drive away from it.

Physical implementation on e-puck robots
The formerly exposed BCM has been implemented on the open-hardware e-puck robot specifically designed for experimental swarm robotics Mondada et al. (2006).The e-puck is a cylindrical robot, of approximately 70 millimeters in diameter, with 8 infrared (IR) proximity sensors (emitter and receptor) located nearly evenly on its perimeter (see Fig. 4).Each sensor is composed of an IR emitter and a receptor, and can be used passively (only the receptor, to sample ambient light), or actively (emitter/receptor, for proximity search or message exchanges).
Our implementation uses only IR sensors (passively and actively, alternately), and the two motors.Each timestep is separated into five phases: 1. Message broadcasting: Each robot first broadcasts a message to all neighboring robots via IR (20-40 cm of range, depending on the angles of the emitter/receptor sensors).The message contains its ID, and the number of agents that have been counted at the previous timestep N ρ (t − 1). 2. Message reading: The received messages are read to measure N ρ (t), the number of detected neighbors.A filtering step also allows detection of robots that were in range in the previous steps, but whose messages were not received (to correct possible miscommunication between robots).A memory of the detection of the other robots in the last 5 steps is maintained, and a majority window algorithm (MWA) is applied to correct for robots not being correctly detected.The filtered value N ρ (t) is thereby constructed.
The MWA counts as a neighbor any robot that has been in range at least 3 times in the last 5 timesteps, or that has just been detected during the present timestep.The number of neighbors can only be underestimated, as we use an unmodulated light source (on a continuous current), and the IR communication between the e-pucks are modulated.This means that no false positives can occur in the communication among the robots, and only false negatives (undetected neighbors) can happen.To try and reduce the influence of non-detection of neighbors, we artificially force the number of neighbors at time t to be bounded below by (N ρ (t − 1) − 1).This mechanism ensures a smooth behavior of the speeds of the motors, by avoiding large jumps in the number of neighbors over time.We define V ρ (t) to be the set of the values N ρ (t − 1) received from the neighbors at time t.
Then the final value N ρ (t) is defined as: 3. Light source sampling: The intensity of the light source is passively sampled on each of the eight IR sensors to obtain S0 (t), S1 (t), ... S7 (t) .The normalized values are then computed as S i (t) = Si (t) / max i Si (t) for i = 0, 1, .., 7. 4. Velocity re-actuation: With the data gathered, the velocity of each motor is computed and updated following Eq.( 12), with α = 100 and β = 20, and a base speed of 100: The speeds here are given in steps/second.1000 steps correspond to one full rotation of a wheel, whose diameter is approximately 41mm.A speed of 500 steps/s corresponds to 500 1000  This allows us to consider the average speeds on the limit cycle as being constants, and to rewrite the BCM on the limit cycle as: A relation between these speeds and the interaxial distance of the robots, leads to the radius of the limit cycle The angular speed can also be expressed as

Dictionary theoretical model-robotic implementation
A formal proof of the convergence of our BCM would be very difficult to provide; this is why we chose to illustrate that our BCM exhibits the same characteristic behavior as the mathematical model of Sect.2, whose convergence we have formally proven.To better link the model and the robotic implementation, Table 1 illustrates the translation of these main characteristics between the theoretical and robotic control mechanisms.The explicit expression of the limit cycle in both cases is also summarized, to show how both depend on their control parameters In the BCM, we see that our initial choice of setting α and β fixes the radius of the limit cycle (and the angular speed on it).If we choose a larger swarm of robots, we should not see a difference in the radius as long as we do not change these two control parameters.Theoretically, we could see the swarm trying to arrange on a circle whose perimeter is too small for all the robots to fit.This means that we should see α and β as functions of N , to be sure that all robots can fit on the consensual limit cycle.We should have the perimeter of the limit cycle at least larger than N times the diameter d of one robot: This result from the BCM matches the mathematical model, where the radius of the limit cycle depends only on the fraction of neighbors M and the radius of observation ρ.In the mathematical model, changing the swarm number does not change the limit cycle's radius either.

Experimental validation
To demonstrate the convergence of our BCM regardless of the initial conditions, tests were carried out with a swarm of eight robots, from 25 different (random) initial situations (positions and headings of the robots).The choice of these initial conditions was made to experimentally validate the inter-and intra-convergence variability.The headings changes in similarly located robots in different initial conditions (e.g., experiments 3 and 17) show how the control mechanism can achieve convergence to a consensual orbit with initial conditions  6, along with the position of the IR source, which remains unchanged during the experiments.Each experiment lasted from 1 to 4 minutes; long enough to let the robots rearrange from their initial positions to the limit cycle, and then maintain the limit cycle.The tracking step allowed us to obtain the position of each robot through all of the 25 experiments.Figure 7 shows the tracking step applied to the videos, while Fig. 8 shows the resulting tracked trajectories.Details on the multi-agent tracking tool developed can be found in Appendix section "Multi-agent tracking tool".

Experiment results
Data extracted from tracking have been used to quantitatively assess convergence, based on the following measures: 1. Mean and standard deviation of each robot's radius in the last 15% of each experiment (empirically found to be in the limit cycle regime for every initial situation).2. Mean and standard deviation of each robot's angular speed in the same time window.The extracted data are shown in Fig. 9 for each experiment, where these averages and standard deviations are pictured for each robot on the left side, and one average for the whole swarm on the It is worthwhile note the very low fluctuations on the average radius and angular speeds for the swarm, through the whole set of initial conditions.
As a measure of the uniformity of the robots' dispersion on the limit cycle, we now recall the order parameter (see Acebrón et al. (2005)) for a group of N oscillators on R 2 in polar coordinates: which takes a real positive value.Note that if all the oscillators share the same phase θ i (t) at a given time, the value for R(t) will be high.Conversely, if all the phases of the oscillators are uniformly distributed over [0 : 2π], R will vanish.
In this paper, we used information on the position in the last 15% of each run to determine the phase of each robot through time, and compute R(t) for the swarm.Figure 10 reports the average in this time window (along with the standard deviation) for each experiment.

Remark:
The positions of robots that got stuck against the IR source during an experiment are shown in the initial configurations plot in Fig. 6, even though their data were not used in the analysis.Over the 25 experiments, 15 robots got stuck, with a maximum of 2 stuck robots for one initial condition (7.5% of unusable data).As the IR source is an attractive point, and because of the avoidance mechanism used, the robots sometimes change their course and end up hitting the light source.A robot in this situation cannot get back to the limit cycle, and stays at the same spot with 0 velocity until the end of the experiment.In the mathematical model of Sect.2, agents can leave the attractor only because of the noise source.Our robots cannot do the same, due to the physical presence of the IR source blocking their path.
This issue could not be solved by using a different avoidance mechanism, as avoiding a light source, solely with the use of the IR sensors is impossible (at short range, the IR sensors are saturated by the light source).This result conforms to our expectations, as it fits the theoretical predictions as well as the underlying Eq. ( 1).

Comparison with the theoretical results
Theoretical results derived in Sect.3.2 lead to the following limit cycle radius and angular speed (with the e-puck robot's interaxial distance e = 52[mm]): These results very closely match our experimental results, where the average radius of the limit cycle over the 25 experiments is L α,β := 252.7 [mm], and the average angular speed is ω α,β := 0.197.

Generalization to other Hamiltonian functions
The choice of a BCM limits the characteristics of the Hamiltonian functions that can be adapted to a hardware implementation.The scope of this paper was not to study the limit of validity of our model; other types of control mechanism may adapt to general Hamiltonian functions.Nevertheless, the BCM described here for the case of the harmonic oscillator (H (x, y) = 1 2 (x 2 + y 2 )) can be extended to more complex Hamiltonian functions.For example, one can obtain Cassini oval-like curves by placing the robots in an environment with two light sources.Using exactly the same control mechanism of Eq. ( 12), we obtain a swarm behavior similar to the mathematical case with the Cassini ovals Hamiltonian function.Depending on the observation range ρ of the robots, cases with one connected consensual orbit or two separated orbits can arise.This time, robots sample the ambient light coming from both light sources, acting as the two focal points from the Cassini Ovals Hamiltonian function.
Tests have not yet been carried out on an actual swarm of robots, but numerical results are shown in Fig. 11.These results are obtained by simulating Eq. ( 12) with dt = 0.1s, on a swarm of 100 robots starting with random positions and headings on the plane.The two light sources are positioned at (±100, 0), and WGN sources are added to the (exact) gathered information to simulate the noise actual robots would experience: Empirically chosen values for the WGN sources are σ l = 0.05 for the light source sampling, and σ N = 0.5 for the number of instantaneous neighbors ( ÑV (t) := N V (t) + σ N dW N (t)).

Conclusion
For a whole class of interacting autonomous agents evolving in an environment subject to noise or external perturbations, we have been able to perform an analytical characterization corroborated by an extensive set of simulation experiments and ultimately to implement the control mechanism on a swarm composed of eight e-puck robots.The possibility to simultaneously complete these complementary approaches on the same class of control mechanisms is exceptional, and originates from the simplicity of the underlying dynamics (i.e., MCD dynamics).However, despite its simplicity, our modeling framework incorporates a number of basic features rendering it generic, namely: non-linearities in the control mechanism; finite range interactions between agents; self-organization mechanisms leading to dynamic consensual spatio-temporal patterns; and environment noise, which corrupts ideal trajectories and sensor measurements.Beyond the pure pedagogic insights offered by this class of models, the results provide strong evidence that built-in stability of the control mechanism and the resilient behavior observed during the actual implementation open the door for more elaborate implementations.The modeling framework exposed here is more than a mere class of models; it offers a constructive method to deduce generalized classes of fully solvable agent dynamics.In particular, extension of the MCD dynamics via the introduction of Hamiltonian functions, involving more degrees of freedom, will allow analysis of the robots' evolution in higher dimension (from planar to 3-D spaces).
Acknowledgments Detailed comments from anonymous referees strongly contributed to the presentation and quality of this paper.Guillaume Sartoretti and Max-Olivier Hongler are partly funded by the Swiss National Science Foundation and the ESF project "Exploring the Physics of Small Devices".Marcelo Elias de Oliveira is supported by the EU-ICT-FET project "ASSISIbf", No. 601074 and by the ESF project "H2Swarm".This last project is funded by the Swiss National Science Foundation under grant number 134317.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Appendix: MCD Dynamics
Consider the dynamics of a harmonic oscillator with Hamiltonian function H (X, Y ) = 1 2 X 2 + Y 2 .The resulting dynamics reads: Clearly, the trajectories solving Eq. ( 20) are circles on the plane R 2 with radius given by the initial conditions (X 0 , Y 0 ).We shall now introduce energy dissipation into the nominally (energy-conserving) Hamiltonian system Eq.( 20).This will be achieved by adding diagonal terms to the matrix V 0 : where U (H ) is a function of the Hamiltonian and we denote U (H ) := d d H U (H ).For instance, the choice U (H ) = −H ⇒ U (H ) = −1 simply yields a linearly damped dynamics.In this case, the trajectories will simply be spirals on R 2 , approaching the origin as time t increases.Observe from Eqs. ( 20) and ( 21) that the conservative part of the dynamical vector field (VF) is orthogonal to the dissipative part: we have the vanishing scalar product (22) The dynamics of Eq. ( 21) will finally be driven by a couple of external noise sources dW X (t), dW Y (t) that will be taken as standard independent WGN processes: Accordingly, the planar dynamics in Eq. ( 23) simultaneously model: (i) a canonical driving, (ii) an energy dissipation mechanism, and (iii) an energy source delivered by the noise sources.The use of WGN implies that X (t) and Y (t) are Markovian processes, for which all probabilistic characteristics of the realizations are entirely given by transition probability densities P(x, y, t | x 0 , y 0 , t 0 ) := P for t > t 0 .Moreover, P(x, y, t | x 0 , y 0 , t 0 ) obeys the FP equation: Asymptotically in time in Eq. ( 23), a balance between dissipation and supply of energy is reached and this leads to a stationary regime, in which all probabilistic features of the dynamics become time invariant.This stationary regime is described by the stationary probability measure P s (x, y) solving the stationary (FP) that results when the left-hand side in Eq. ( 24) equals zero.The cylindrical symmetry of Eq. ( 23) implies that, in the stationary regime, all probabilistic features of the model must themselves be cylindrically invariant; hence, P s (x, y) = P s [H (x, y)] = P s (H ).By direct calculation, we can indeed verify that (with N the normalization factor).( 25) Now, provided that Eq. ( 25) is integrable, meaning that N < ∞ (i.e., implying an overall dissipative dynamics), the stationary probability characteristics of the dynamics are given by Eq. ( 25).In particular, for the linearly dissipative case obtained when U (H ) = −H , we consistently obtain a bivariate Gaussian density with cylindrical symmetry: .

Multi-agent tracking tool
To determine the states of every single agent in our experiments, we developed a multiagent tracking tool based on computer vision techniques and the Kalman filter.This task was particularly challenging due to the very inhomogeneous lighting conditions in the acquired videos (mainly because of the central light source).This acquisition setup is composed of an overhead digital camera equipped with a wide angle lens (Canon EF-S 10-22mm f/3.5-4.5 USM, mounted on a Canon EOS 7D body).Our multi-agent tracking tool has been developed in C++ using the GNU compiler collection (GCC 4.6.3)and standard computer vision algorithms included in the open source computer vision library (OpenCV, v. 2.3.1-7,http://opencv.org/).The employed technique follows three steps: 1.The image foreground is extracted for further processing using a precomputed mean background.2. The Gaussian scale-space representation of the extracted foreground I F (x, y) containing the agents was convolved with a Gaussian kernel, as L(x, y, σ ) = G(x, y, σ ) * I F (x, y), where, in a two-dimensional space, the Gaussian function is given by G(x, y; σ ) = 1 2πσ 2 e −(x 2 +y 2 )/2σ 2 (26) and σ is the standard deviation.The outcome of this convolution at a specific scale is strongly dependent on the relationship between the sizes of the blob structures in the image domain.Since our background is characterized by high pixel values and the blobs dimensions are known, every agent can be characterized by a local extrema of the resulting Gaussian scale-space representation L(x, y, σ ). 3. Assuming that the initial frame does not present any false positives, every single detected blob in the initial frame represents an agent (see Fig. 7, left).The detected spatial locations of the agents are used as inputs for the Kalman filter Kalman (1960), where the initial speeds of the agents are assumed to be zero.The estimates of the system states are based on the initial conditions and the sequence of the detected blobs along the experiments.Since the radius of the cylindrical robots are known, as well as their directions and linear speeds, a new measurement is easily casted as a good candidate input for the Kalman filter or as a new input for the designed regressive filter.
In this work, the agents are assumed to follow non-uniform translational motions, with a constant acceleration between consecutive frames.The uncertainties associated with the measurements of the positions induced by different lighting conditions are treated as white noise.The states x i k of the i-th agent at time k as well as the measurement z i k are described by linear SDE in the image domain: where x i k ∈ 4 , z i k ∈ 2 .w k−1 and v k are WGN sources, respectively, representing the process and measurement noises.They are assumed to be mutually statistically independent, and to be drawn from zero-mean multivariate Gaussian distributions with covariances Q and R. The covariances' matrices can be described as a function of the experimental parameters, such as the lighting variability across the experiments, the number of agents, etc.However, in our experiments, they were assumed to be constant.The evolution of the current state vector x i k over time is described by the transition matrix derived from the non-uniform translational motion equations where t is equal to 0.1 s, as the videos were acquired at a constant frequency equal to 10 Hz, the time elapsed between consecutive frames was assumed to be constant.H is a 2 × 4 transformation matrix that related the predicted state vectors x i k into the measurement domain: Remark: The videos of all the experiments can be found online, at http://www.youtube.com/channel/UCCgwiDNHoyG9RJ2_pdyewmg.

Fig. 2
Fig. 2 Numerical simulation for N = 200 agents, ρ = 1, σ = 0.07 and M = 4.The black dots show the final position of the agent at time T = 50[s], while the blue line gives the trajectory of one arbitrary agent (with its observation range in red).The black circle (hidden below the swarm of agents) gives the theoretically computed limit cycle.This shows the exactitude of Eq. (10), even for finite swarms.Left Harmonic

Fig. 4
Fig. 4 Left Theoretical robot, seen from above.Note the eight sensors evenly distributed on the perimeter of the robot's body and the two motors L and R. The point S (resp.U ) shows the stable (resp.unstable) facing direction of the light source in the BCM.Right Actual e-puck robot in the same position; the eight proximity sensors (ps0,ps1,...,ps7) are nearly evenly distributed on the perimeter

Fig. 5
Fig. 5 Organigram of the robot's main programming loop

Fig. 6
Fig. 6 Initial positions and headings of the eight robots in each of the 25 experiments.Each circle corresponds to a robot, and the corresponding line denotes the heading.The center cross represents the position of the light source

Fig. 7 Fig. 8
Fig. 7 Left A randomly selected initial condition assuming eight aligned equally spaced agents in the inferior working space region.The IR radiation source is located at the center of the rectangular arena.Right Respective agents' control mechanism and their mutual interactions.It is important to observe that the proposed Kalman filter is capable of determining the position of an agent with high precision independently of lighting conditions and other agents' positions

Fig. 9
Fig. 9 Average radius/angular speed in the last 15 % of each experiment.Left For each robot.Right Averaged for the whole swarm

Fig. 11
Fig. 11Numerical simulation of the BCM on a swarm of 100 robots, with two light sources located at (±100, 0).With an observation range of ρ = 100, the swarm consensually rotates around both light sources (left), whereas a smaller range of ρ = 10 lets the swarm split into two sub-swarms rotating around each light source (right).The trajectory of 10 randomly selected robots are shown, along with the final position of all robots as black dots 2 are the positions of the planar agent, H (X, Y ) : R 2 → R + and U [H ] : R + → R (U [H ] stands for the derivative of U [H ] with respect to H ) are differentiable functions, σ > 0 and ω ∈ R are constant parameters.The inhomogeneous terms dW X (t) • π • 41 ≈ 64.4 [mm/s].5. Avoidance mechanism: A proximity search is performed on each sensor to detect obsta- cles; avoidance is implemented by adapting the motors' velocities.The avoidance mechanism used is a standard Braitenberg-inspired avoidance mechanism.It has been designed to avoid collisions between e-puck robots, but cannot guarantee collision avoidance with the light source.Its pseudo-code is given in Algorithm 1.

Table 1
Theoretical and robotics implementations of the main control mechanisms in our model, and their respective control parameter (s)