Stigmergy: from mathematical modelling to control — TEI-tagged view

Alain BoldinipersName, Martina CivitellapersName, Maurizio PorfiripersName (2024date). “Stigmergy: from mathematical modelling to control.” Royal Society Open Sciencetitle·j.

Object of study — a continuum model linking the density of environmental tracesterm·dom to the density of a swarmterm·key, so that stigmergyterm·dom can be designed

Research question — which modifications of the environment are needed for a swarmterm·key to reach a desired behaviour

Contribution — a continuified view of stigmergyterm·dom that computes explicitly the traceterm·dom distribution needed for a desired swarmterm·key configuration

Full text — tagged

Stigmergyterm·dom: from mathematical modelling to control Abstract Stigmergyterm·dom, the indirect communication between agents of a swarmterm·key through dynamic environmental modifications, is a fundamental self-organization mechanismterm·dom of animal swarmsterm·key. Engineers have drawn inspiration from stigmergyterm·dom to establish strategies for the coordination of swarmsterm·key of robots and of mixed societies of robots and animals. Currently, all models of stigmergyterm·dom are algorithmic, in the form of behavioural rules implemented at an individual level. A critical challenge for the understanding of stigmergic behaviour and translation of stigmergyterm·dom to engineering is the lack of a holistic approach to determine which modifications of the environment are necessary to achieve desired behaviours for the swarmterm·key. Here, we propose a mathematical framework that rigorously describes the relationship between environmental modifications and swarmterm·key behaviour. Building on recent strides in continuificationterm·key techniques, we model the swarmterm·key and environmental modifications as continua. This approach allows us to design the environmental modifications required for the swarmterm·key to behave as desired. Through analytical derivations and numerical simulations of one- and two-dimensional examples, we show that our framework yields the distribution of tracesterm·dom required to achieve a desired formation. Such an approach provides an adaptable framework for different implementation platforms, from robotic swarmsterm·key to mixed societies of robots and animals.

Keywords: animal behaviour, collective behaviour, control, swarmsterm·key Introduction Stigmergyterm·dom—defined as a ‘mechanism of indirect coordination in which the traceterm·dom left by an action in a medium stimulates subsequent actions’ [ 1 ]—is a fundamental means of self-organization in complex systems [ 2 ]. The concept of stigmergyterm·dom originated from the study of coordination in groups of animals [ 3 ]. For example, ants release pheromones in the environment [ 4 ] to guide other ants back to the nest once they find food [ 5 ], or to cooperatively transport large food items [ 6 ]. Likewise, pumas and other large Felidae leave tracesterm·dom in the environment to mark territory so that other conspecifics would avoid the area [ 7 ]. In animal colonies that build structures [ 8 , 9 ], such as wasps [ 10 ] and termites [ 11 , 12 ], stigmergyterm·dom allows the incredible coordination of thousands of individuals to build intricate nests, without any pre-planning and central coordination.

With the advent of robotic systems, engineers have taken inspiration from stigmergyterm·dom in natural swarmsterm·key to build self-organizing groups of robots that coordinate to achieve a desired goal [ 13 – 15 ]. In stigmergyterm·dom, robots in a swarmterm·key not only coordinate to sense the environment [ 16 , 17 ], but they actively modify it to communicate with each other. Through stigmergyterm·dom, robotic swarmsterm·key can gather objects [ 18 ], sort them [ 15 ], navigate unknown environments [ 19 ], and search and track a moving target [ 20 ]. Stigmergyterm·dom is particularly promising for collective constructionterm·key, as the structure that is being built can be utilized as the stigmergic signal itself [ 21 , 22 ]. In this vein, extensions of stigmergic approaches have been proposed to achieve precise and accurate collective constructionterm·key [ 23 , 24 ]. Stigmergyterm·dom is an important coordination strategy not only in biological and robotic systems but also in so-called ‘mixed societies’ that integrate biological and robotic individuals [ 25 , 26 ]. Therein, robots are used to create new collective responses or to elicit a desired behaviour from the animal group [ 27 ]. In addition to direct social interactions between animals and robots, stigmergic stimuli released by robots have been often used to achieve the desired swarmterm·key behaviour, for example through light signals [ 28 ] or pheromone release [ 29 ] by robots.

A fundamental challenge for biologists is to reconstruct how environmental modifications implemented at an individual level, without a central coordination, trigger a complex behaviour of the swarmterm·key, such as the construction of a nest or food gathering. The literature relies on agent-based modelling grounded in individual behavioural rules [ 13 – 15 ], which do not allow a holistic analysis of the swarmterm·key. Currently, no stigmergic model allows the study of the overall distribution of environmental modifications that enable these complex behaviours. Similarly, engineers are faced with the design problem of how to control robotic agents toward achieving a desired behaviour for the swarmterm·key, be it robotic or mixed. Such a challenge requires the development of a mathematically tractable and interpretable framework for stigmergyterm·dom that could complement existing algorithmic implementations, which rely on a sequence of discrete behavioural rules implemented at the individual level.

Here, we put forward a control-oriented mathematical backdrop to holistically describe environmental modifications in animal swarmsterm·key and engineer stigmergic interactions in robotic and mixed swarmsterm·key, where units of the swarmterm·key interact with each other and with dynamic modifications of the environment that we identify as ‘tracesterm·dom’. Our mathematical framework is based on recent advancements in the field of control of large multi-agent systems, which relies on a ‘continuificationterm·key’ of the equations of motion of the swarmterm·key [ 30 , 31 ] (see figure 1 ). In this vein, rather than modelling and controlling individual motions [ 32 ], we model the swarmterm·key as a continuous fluid and describe the spatio-temporal evolution of its density, similar to thermodynamic approaches for large-scale systems [ 33 ] and former robot density control algorithms [ 34 – 38 ]. Our work differs from previous studies on continuification-based control of large swarmsterm·key [ 30 , 31 ] as we do not apply a control action to each swarmterm·key unit. In animal swarmsterm·key and mixed societies, it is not possible to directly control biological units, and even in robotic swarmsterm·key controlling each unit would require a centralized approach or, at least, some form of connectivity in the swarmterm·key [ 39 ]. Based on these groundings, we focus on the realistic case where the control engineer can only design dynamic modifications of the environment, such that no centralization or communication between the units is needed.

Schematics of the control problem and proposed framework for our mathematical model of stigmergyterm·dom. The control problem that we seek to address is achieving a desired configuration of a large swarmterm·key in terms of position of its discrete units. Such a desired configuration is continuified (1) to provide the desired density of the swarmterm·key. The density of the tracesterm·dom is computed (2) based on the desired density of the swarmterm·key, and is then discretized (3) to allow actual deployment of the tracesterm·dom. The final, desired configuration of the swarmterm·key arises from the interaction between units of the swarmterm·key and the tracesterm·dom left (4).

While our work shares the same backdrop of continuified swarmsterm·key as that in Maffettone et al . [ 30 ], it tackles the control problem through a totally different perspective. Our control action is not directly applied to the swarmterm·key agents; rather, it is mediated by the deposition of tracesterm·dom in the environment. Such a continuified view of stigmergyterm·dom is an original contribution of this paper, which allows us for the first time to explicitly compute the traceterm·dom distribution necessary to achieve a desired swarmterm·key configuration. In this vein, we offer closed-form solutions of the control input in the form of density of tracesterm·dom.

We examine a prototypical example of a task that can be achieved through stigmergyterm·dom: shape formation [ 21 , 40 ]. We seek to compute the distribution of tracesterm·dom that make the swarmterm·key achieve a desired formation, that is, a prescribed density profile. To this end, we propose a general procedure to compute the density of tracesterm·dom that allows us to retrieve the prescribed density of the swarmterm·key. Once the density of tracesterm·dom is known, it is discretized to allow for deployment in a real environment, in a discrete form. We do not focus on implementation details of how such density profile of tracesterm·dom can be generated, as we concentrate on the mathematical modelling of the overarching stigmergic control strategy, generalizable to a variety of practical implementations. While one could describe the dynamics of deposition of the tracesterm·dom, such a process would be heavily context-dependent (for example, the dynamics of ants’ pheromone release [ 13 ] is different from that of robotic collective additive manufacturingterm·dom (AMterm·dom) [ 23 ]), thus hindering the generality of our theory.

We demonstrate the potential of our approach through a series of one- (1D) and two-dimensional (2D) case studies. In the 1D benchmarks, we show that our approach allows for the swarmterm·key to achieve a stationary formation and generate travelling wavesterm·key [ 41 ]. The 2D simulations validate our approach in a more realistic scenario. We select a historical, fascinating example, by replicating the complex, non-uniform motion of the robotic lion of Leonardo da Vinci along a wavy circle [ 42 ] with our swarmterm·key.

Mathematical model We start by introducing the general framework of our problem in a 1D setting, considering a swarmterm·key coordinating over a ring (§2.1). We define two types of interactions between units on the ring (§2.2): a purely repulsive interaction and an attractive–repulsive interaction. In §2.3, we put forward analytical solutions for both of these types of interactions. We then introduce a numerical method to study the time-evolution of the swarmterm·key at the continuum level in the 1D problem (§2.4), along with the discrete counterpart used to validate our approach (§2.5). Finally, in §2.6, we extend our mathematical model to multi-dimensional problems.

General framework in one dimension Our aim is that of controlling a swarmterm·key moving on a ring ( S 1 ) [ 43 ], parametrized with the angle x ∈ [ − π , π ) . The state of the swarmterm·key is described by the density ρ : 𝕊 1 × ℝ + → ℝ + . We assume that the evolution in time of the density is governed by mass conservation: where ( ⋅ ) t and ( ⋅ ) x indicate partial derivatives in time and space. Here, v : 𝕊 1 × ℝ + → ℝ is the velocity of the swarmterm·key.

We adopt a kinematic model [ 44 ] that neglects inertia effects, such that the velocity v depends explicitly on the density through a circular convolution [ 45 ]: where q T : 𝕊 1 → ℝ is an interaction kernel and f : 𝕊 1 × ℝ + → ℝ is an external forcing term. The circular convolution is equivalent to the linear convolution over ℝ of a periodic function defined over 𝕊 1 with a non-periodic kernel defined over ℝ (see §2.2). The kernel embodies the interactions between swarmterm·key units, typically derived from a potential, which may be repulsive, attractive or a combination of the two modulated by the distance [ 43 ]. By construction, the circular convolution in equation (2.2) guarantees that lim x → π ⁡ v ( x , t ) = v ( − π , t ) , so that ρ ( − π , t ) v ( − π , t ) = lim x → π ⁡ ρ ( x , t ) v ( x , t ) . Thus, the total mass of the swarmterm·key ∫ − π π ρ ( y , t ) d y = M is conserved. The model in equations (2.1) and (2.2) naturally arises from the continuificationterm·key of classical, discrete kinematic models that have been widely used to describe agents’ behaviour in swarmsterm·key [ 43 ].

We seek to develop a control strategy where the external forcing is provided by some ‘tracesterm·dom’ with density ρ tr : 𝕊 1 × ℝ + → ℝ , interacting with the swarmterm·key, such that Since f is defined as a circular convolution, lim x → π ⁡ f ( x , t ) = f ( − π , t ) by design. Here and in the following, we adopt the same interaction kernel q T for the interaction between units in the swarmterm·key and the tracesterm·dom. Should the units of the swarmterm·key be able to distinguish between their peers and the tracesterm·dom, different interaction kernels can be utilized. We accept the possibility that ρ tr may be negative as a computational device to simplify our derivations; we can recover a non-negative ρ tr with equivalent effects, as made clear in §2.3.

Periodicization of the interaction kernels Periodic interaction kernels q T are obtained from the periodicization of standard non-periodic kernels q : We consider the following two standard non-periodic kernels.

Repulsive kernel. The non-periodic repulsive kernel is in the form where sgn is the sign function. Note that we utilize a length-scale L while fixing the domain to [ − π , π ) . Periodicization leads to By separating the infinite series into two other infinite series based on the sign of x + 2 k π and computing each of these series individually leads to Morse kernel. The non-periodic Morse kernel is an attractive–repulsive kernel, in the form having defined By following similar steps to the repulsive kernel, we find where Solution in one dimension Let us assume that our objective is to attain a chosen density ρ ‾ ( x , t ) : 𝕊 1 × ℝ + → ℝ + for the swarmterm·key. To find the forcing term f necessary to obtain the desired density, we assume that ρ ‾ satisfies equation (2.1) , such that with v ¯ ρ ( x , t ) = ∫ − π π q T ( x − y ) ρ ¯ ( y , t ) d y . By integrating equation (2.12) , assuming ρ ‾ ≠ 0 for any x , we obtain with A ( t ) being a function of time only. We note that A ( t ) may not be chosen freely, as f should satisfy the integral condition ∫ − π π f ( x , t ) d x = ∫ − π π q T ( x ) d x ∫ − π π ρ t r ( x , t ) d x , according to Fubini’s theorem [ 46 ]. As typically we consider odd kernels q T , we must have that f integrates to zero over the domain. Thus, A ( t ) is chosen to instantaneously satisfy this condition.

Once f is computed in closed form, we can derive ρ tr from the deconvolution of equation (2.3) . Such an operation depends on the interaction kernel q T . For standard interactions kernels introduced in §2.2, we have established a closed-form analytical solution.

Repulsive kernel. For the repulsive kernel [ 47 ], we can expand the convolution in equation (2.3) as By taking two derivatives with respect to x , we find Thus, we can retrieve by integration where B ( t ) is an arbitrary function of time.

Morse kernel. With the Morse kernel, we can write the convolution in equation (2.3) as where From equation (2.16) , we find However, from equation (2.13) , we only have access to f , not to f R and f A . From differential and algebraic manipulations of equations (2.17) and (2.19) , we find a problem in f R and f only: Equation (2.20) can finally be solved for f R through the use of spatial Fourier series.

The distribution of tracesterm·dom that generates the desired swarmterm·key density is not unique. In fact, the solution is defined up to an additive term ρ ~ tr ( x , t ) such that, for any x , In particular, any static equilibrium ρ eq ( x ) of equation (2.1) , with corresponding velocity v e q ( x ) = ∫ − π π q T ( x − y ) ρ e q ( y ) d y = 0 for any x , is a solution of equation (2.21) . For example, for a purely repulsive kernel with interaction length L and a static desired density, we obtain where B is an arbitrary constant.

Numerical method To study the time-evolution of the equations of motion, we derive a numerical solution. Specifically, to solve the hyperbolic equation in equation (2.1) , we rely on a finite volume methodterm·key, which naturally accounts for conservation of a given quantity and for the finite velocity of perturbations [ 48 ]. We divide the domain [ − π , π ) in N C cells C i ( i = 1 , … , N C ) with uniform size Δ x and consider as variables the average value of the density over the cells at time t n : The evolution of this quantity in time is governed by where Δ t is the integration time-step and F i + 1 2 n and F i − 1 2 n are the numerical fluxes on the right and left interfaces, which approximate the value of the flux f ( ρ ) = ρ v ( ρ ) at the interact over the time-step: For simplicity, we consider the Lax–Friedrichsterm·key method, which reads In 1D simulations, we utilize N C = 500 cells and a time-step of Δ t = 0.1 . To quantify the error between actual density of the swarmterm·key and desired one, we leverage the Kullback–Leiblerterm·mth (KLterm·mth) divergence of the actual density of the swarmterm·key from the desired one over time, which represents a form of statistical distance [ 49 ].

Discrete framework All the results in the previous subsections are only valid at the continuum-level description of the swarmterm·key. To validate our overall approach in more realistic conditions, we conduct 1D simulations of the discrete swarmterm·key. We consider the problem of coordinating N units of a swarmterm·key, each with mass m = M / N . The dynamics of the units is described by their position x i ( t ) along the circle, and is governed by a discretized version of equations (2.1) , (2.2) , and (2.3) , namely [ 43 ], with f i being the external forcing on x i . This term is provided by the interaction of the units with the N tr discrete tracesterm·dom at positions x i tr : The desired configuration of the swarmterm·key is imposed in terms of density, so that an initial continuificationterm·key is not required. Once the desired density of the tracesterm·dom is computed, self-equilibrating solutions are added to ensure that the density is positive everywhere. The resulting distribution of tracesterm·dom is discretized over a certain, fixed number of cells, such that tracesterm·dom can be added to the simulation. The number N tr of tracesterm·dom is a result of the algorithm, as we approximate it as the closest integer to ( ∫ − π π ρ t r d x ) / m . Upon simulating the dynamics of the discrete system, we can define the corresponding continuum counterpart through continuificationterm·key, which is achieved from kernel density estimation with a manually tuned Gaussian estimation kernel [ 50 ]. In this case, the KLterm·mth divergence is evaluated on the continuified distribution.

Extension to multi-dimensional problems Our method can be readily extended to multi-dimensional problems. Specifically, we consider a density ρ : 𝕊 n × ℝ + → ℝ + , where n = 2 or n = 3 and 𝕊 n is the n -sphere parametrized with [ − π , π ) n , and extend equation (2.1) to with 𝐯 : 𝕊 n × ℝ + → ℝ n being the velocity vector and ∇ ⋅ ( ⋅ ) the n -dimensional divergence operator. Equation (2.2) is generalized as a multi-dimensional circular convolution with vectorial interaction kernel 𝐪 T : 𝕊 n → ℝ n and external forcing term 𝐟 : 𝕊 n × ℝ + → ℝ n : The formula for 𝐟 is a multi-dimensional version of the circular convolution in equation (2.3) , with a density of tracesterm·dom ρ tr : 𝕊 n × ℝ + → ℝ .

In the multi-dimensional case, closed-form solutions to the problem are not available. Specifically, closed-form solutions to the periodicization of the interaction kernels in all spatial directions cannot be easily found. As such, we truncate the multi-dimensional infinite series extending equation (2.4) to a finite number of terms. For example, for n = 2 , where 𝐪 ( x , y ) is a non-periodic interaction kernel in 2D.

To obtain ρ tr , we numerically solve an inverse problem for equation (2.29) , substituting the density ρ of the swarmterm·key with the desired density ρ ‾ . In our 2D simulations, we utilize a finite volume approximation of the equation with a complete 2D Lax–Friedrichsterm·key flux. We utilize 30 cells to discretize each direction, for a total of 900 cells, along with a time-step of Δ t = 0.1 .

Results We first demonstrate our approach on 1D problems over a ring, for which we propose a detailed analysis of the density of tracesterm·dom required to achieve a desired density of the swarmterm·key, along with the error over time (§3.1). Further, we validate our approach by simulating the same 1D problems with discrete swarmsterm·key and tracesterm·dom. Then, in §3.2, we consider a 2D example of a complex, evolving formation problem of a swarmterm·key that can be achieved with our stigmergic approach.

Results in one dimension We consider three applications of our approach in 1D: (i) a swarmterm·key that reaches a stationary formation; (ii) a swarmterm·key that forms a travelling waveterm·key in a non-dispersive medium, where the group and phase velocities are equal; and (iii) a swarmterm·key that forms a travelling waveterm·key in a dispersive medium, where the group and phase velocities are different.

As the waveform, we consider a von Mises function with mean μ and concentration coefficient k [ 51 ]: with J 0 being the modified Bessel function of the first kind of order zero. In the case of a travelling waveterm·key in a non-dispersive medium, we let the function travel at a constant speed v , without affecting the shape of the wave, such that the mean varies with time as μ + v t . Finally, for the simulation of a travelling waveterm·key in a dispersive medium, we let k ( t ) = k 0 + A sin ⁡ ( ω t ) to capture oscillations around a fixed value k 0 , with amplitude A and radian frequency ω .

We run a simulation for each of the three applications, utilizing a finite volume methodterm·key to evolve the density of the swarmterm·key over time, as described in §2.4. All simulations start from a uniform density distribution with unitary mass. For all simulations, we consider a repulsive kernel with interaction length L = 1 , set μ = − 1 , and run the simulation until t F = 100 . For the stationary formation and the travelling waveterm·key in a non-dispersive medium, we fix k = 2 . In both travelling waveterm·key simulations, we set v = 0.05 . For the case of the travelling waveterm·key in a dispersive medium, we choose k 0 = 2 , ω = 0.1 and A = 0.5 , with an initial phase of − π . Results of these simulations are shown in figure 2 , where we also consider the absolute error between actual density of the swarmterm·key and desired one and the KLterm·mth divergence of the actual density of the swarmterm·key from the desired one over time. Videos of the simulations are available as electronic supplementary material.

Results of one-dimensional, continuum-level simulations. ( a,c,e ) The profiles of the density of the swarmterm·key at the end of the simulation ( t = 100 , blue solid line), along with the instantaneous desired profile of the density of the swarmterm·key (red dashed line), for the cases of stationary formation ( a ), travelling waveterm·key in a non-dispersive medium ( c ), and travelling waveterm·key in a dispersive medium ( e ); insets indicate the corresponding density of the tracesterm·dom. ( b,d,f ) The absolute error at each point in the spatial domain over time, for the cases of stationary formation ( b ), travelling waveterm·key in a non-dispersive medium ( d ), and travelling waveterm·key in a dispersive medium ( f ); insets indicate the KLterm·mth divergence of the actual density of the swarmterm·key from the desired one over time.

In the case of stationary formation ( figure 2 a ), the swarmterm·key is able to almost perfectly achieve the desired density ( K L < 0.003 at the end of the simulation). As one may intuitively guess for a repulsive potential, tracesterm·dom have a higher density where the density of the swarmterm·key ought to be small. In particular, we find that the profile of the density of the tracesterm·dom closely resembles the opposite of the desired density of the swarmterm·key, see equation (2.16) . Notably, we observe that there are areas where the density of tracesterm·dom is negative. The indeterminacy in the distribution of tracesterm·dom allows us to obtain a non-negative density. As clear from equation (2.16) , the actual density of the tracesterm·dom can always be made positive from the addition of a constant (that is, a self-equilibrating field of tracesterm·dom). From the analysis of the error over time ( figure 2 b ), we find that the error at the peak quickly decays to zero (more than 90 % of its initial value within 10 time units), while the error in the tail of the function has a longer decay time (with similar decays in 50 time units).

Also when considering the travelling waveterm·key in a non-dispersive medium ( figure 2 c ), we observe an almost perfect tracking of the desired formation ( K L < 10 − 5 after 60 time units). In this case, the density of the tracesterm·dom does not resemble the opposite profile as the instantaneous desired density, as it is corrected to anticipate the change in the desired density over time. Interestingly, we find that the KLterm·mth divergence decays slightly faster than for the stationary formation, reaching K L < 0.1 in 12.4 time units against 14.4 time units for the static case ( figure 2 d ). This result may be explained by the faster decay in the error of the tail that follows the peak (initially on its left). This analysis is corroborated by the error over time, wherein one of the tails decays faster to zero error ( 90 % decay from its initial value in about 20 time units) compared with the stationary formation in figure 2 b , while the other translates at a velocity dictated by v .

For the travelling waveterm·key in a dispersive medium ( figure 2 e ), we record similar profiles for the density of the swarmterm·key and tracesterm·dom as for the non-dispersive medium. The main difference is in the width of the peak of the distribution, which is wider for the dispersive case due to the modulation and the initial different phase. Despite the additional modulation of the waveform, we still obtain similar profiles for the error ( figure 2 f ).

In addition to these simulations at a continuum level, we conduct analogous analyses with discrete swarmsterm·key and tracesterm·dom (see §2.5). In particular, we replicate the three 1D studies at a discrete level. We consider very large swarmsterm·key, with 5000 units, which would be complex to control with traditional strategies [ 52 ]. In figure 3 , we show the KLterm·mth divergence for the discrete simulations, computed from the continuified density of the discrete swarmterm·key: the three discrete swarmsterm·key can all track their desired distribution. These simulations confirm the applicability of our approach in real, discrete swarmsterm·key.

Results of 1D, discrete-level simulations. Blue solid, red dashed and black dash-dotted lines show the KLterm·mth divergence of the actual density of the swarmterm·key from the desired one over time for the three cases of stationary formation, travelling waveterm·key in a non-dispersive medium, and travelling waveterm·key in a dispersive medium, respectively.

Results in two dimensions Our approach can be easily extended to 2D and three-dimensional (3D) problems, as described in §2.6. Herein, we consider 2D simulations to show realistic applications of our method. Specifically, we consider the problem of a swarmterm·key inspired by Leonardo da Vinci’s lion robot [ 42 ]. This ‘automaton’ was designed by Leonardo to autonomously move on a trajectory, stop in front of Francis I King of France, opening its breast full of flowers, and return to its initial position. We use these steps as inspiration for a model example of coordination in 2D, where a swarmterm·key tracks and deforms over a wavy circle, similar to a flower with six lobes; see electronic supplementary material and figure 4 . Specifically, we task the swarmterm·key to aggregate in a compact formation ( figure 4 , bottom) and reach a target location by moving along the wavy circle ( figure 4 , right). Upon reaching the target location, the swarmterm·key shall change the formation (simulating the lion opening its breast; figure 4 , top) to perform a collective task (for example, collective constructionterm·key). Upon completion, the swarmterm·key shall resume its original formation (representing the lion closing its breast) and return to the starting location by moving on the wavy circle ( figure 4 , left). Videos of the simulation are in the electronic supplementary material.

Results of 2D, continuum-level simulations. The figure shows a few snapshots of the density of the swarmterm·key, after a few time-steps ( t = 1 , bottom), halfway through the first half of the wavy circle (right), as the breast of the lion is completely open (top), and halfway through the second half of the wavy circle (left). We show both 3D and top views. In the top view, the red circle indicates the desired trajectory of the centre of the swarmterm·key. The plot at the centre represents the KLterm·mth divergence of the actual density of the swarmterm·key from the desired one over time.

In these simulations, for simplicity, we assume quasi-static changes for the desired density, whereby we assume ρ ‾ t = 0 in our computations. Similar to the 1D case, we start from a uniform distribution with unitary integral. We set the final time to t F = 300 , and divide this time-span in three phases with equal duration: translation along the first half of the wavy circle ( 0 − 100 ), opening/closing of the breast ( 100 − 200 ) and return to the initial position along the second half of the wavy circle ( 200 − 300 ). We require the swarmterm·key to rigidly translate along the trajectory while maintaining a formation described by a 2D von Mises function with k x = k y = 2 . To design the desired trajectory, we let the mean of the 2D von Mises function translate along the wavy circle. The opening and closing of the lion’s breast is simulated by dynamically varying k x of the 2D von Mises function from 2 to 3 .

Our results indicate that our approach is successful in replicating the motion of Leonardo’s automaton ( figure 4 ). Based on the KLterm·mth divergence ( figure 4 , centre), the density of the swarmterm·key quickly reaches the desired one, with overall small differences ( K L < 0.3 after 11 time units). After the initial settling, the KLterm·mth divergence is larger in the first and third phases (around 0.25 ), while it decreases in the second phase during the opening and closing of the breast (around 0.10 ). In principle, adopting a non-quasi-static approach can decrease the KLterm·mth divergence, paying the price of a higher computational cost.

Discussion and conclusions Complex systems composed of many units often display surprising self-organization phenomena [ 53 ]. Self-organization is inevitably associated with information transfer between the units of the system, which occurs through direct social interactions or indirect cues mediated by the environment [ 54 ]. Stigmergyterm·dom encapsulates the ability of the units of a system of transferring information among them by modifying the surrounding environment [ 1 ]. This concept has found a broad range of application in animal [ 2 ] and robotic [ 18 ] swarmsterm·key, mixed societies of animals and robots [ 26 ], and even human social systems [ 55 ].

A particularly challenging endeavour in biology is the analysis of how stigmergic signals left in the environment lead to self-organization and execution of complex tasks, even in the absence of central coordination. Such a challenge persists in engineering design of stigmergic modifications of the environment in order for the swarmterm·key to behave as desired. This problem requires the formulation of mathematically tractable and interpretable models to be used in the analysis of swarmsterm·key and the design of control systems. However, literature on stigmergyterm·dom typically focuses on agent-based models, where each unit of the swarmterm·key follows a set of behavioural rules. These models are difficult to use for holistic analyses of environmental modifications and control purposes.

Here, we propose a new mathematical framework to study stigmergyterm·dom in swarmsterm·key. Through this framework, we lay out how tracesterm·dom should be left in the environment to enable the coordination of a swarmterm·key, provided a mathematical model for the motion of the swarmterm·key is available. Such a framework is rooted in recent ideas for the control of large multi-agent systems [ 30 , 31 ], which draw on the analogy between robots and particles in a fluid. While for small groups of particles we can describe the motion of each individual unit, when considering large systems a continuum description is preferable. Thus, we continuify the problem by modelling the swarmterm·key as a fluid, design the control at a continuum level in terms of density of tracesterm·dom, and discretize this distribution to deploy them in the environment. This approach transforms the control problem on a large number of coupled ordinary differential equations to that on a single partial differential equation, which becomes analytically tractable.

One of the potential applications of the proposed setting is collective constructionterm·key in robotic swarmsterm·key through collective AMterm·dom, where buildings and infrastructure are built through layer-by-layer deposition [ 56 – 59 ] by a team of robots [ 60 ]. This paradigm has been inspired by other efforts in collective constructionterm·key with teams of robots using pre-made parts [ 61 , 62 ], amorphous deposition [ 63 ], or the robots themselves [ 64 ]. While initial efforts focused on ground-based robots [ 65 – 67 ], the recent breakthrough in Zhang et al. [ 68 ] demonstrated the possibility of aerial AMterm·dom for large-scale buildings and infrastructure, going far beyond previous efforts that only focused on structure assembly through drones [ 69 ]. Robots in construction teams ought to reach a specific location and coordinate with each other. Our work provides a pathway to achieve these goals, making a first step toward collective AMterm·dom algorithms.

Within the context of collective AMterm·dom, tracesterm·dom may be left in the environment by the swarmterm·key or by one or more exogenous robots, not directly involved in the construction. While stigmergic cues could be provided by the built structure itself [ 22 , 23 ], it is tenable that in many applications the traceterm·dom deposition process should be independent from the printing. A potential implementation of the traceterm·dom deposition entails the use of a separate nozzle, utilizing a coloured dissoluble material that can be identified by other robots. Tracesterm·dom can be left as signs on printed structures or as thin marks on the ground and removed after construction, similar to supporting material used in desktop 3D printers [ 70 ]. Beyond collective AMterm·dom, the tracesterm·dom needed for our algorithm may assume different forms. For example, they may represent actual fixed, physical objects, an additional swarmterm·key of mobile robots, spatially extended continuous robots that can actively deform or other environmental modifications that affect the swarmterm·key (such as changes in background lighting in light-foraging swarmsterm·key [ 71 ]).

Regardless of the application, we assume that units in the swarmterm·key are not able to distinguish between other units and tracesterm·dom. Such a possibility is in practice verified in large swarmsterm·key where individual units have limited computing power and may only perform simple operations to avoid collisions with surrounding entities, be they other units or tracesterm·dom. When a unit can discern tracesterm·dom from other units, we have an additional degree of freedom, as we can modify the type of interaction units have with tracesterm·dom. For example, should units be able to distinguish different types of tracesterm·dom, negative densities of tracesterm·dom may be implemented in practice by using tracesterm·dom encoding interactions via an attraction potential.

Interestingly, we find that the distribution of tracesterm·dom has several degrees of indeterminacy. Such indeterminacy is not a drawback, but rather an advantage of the method, as it provides additional degrees of freedom for the control and swarmterm·key designer. For example, if tracesterm·dom represent fixed, physical objects that cannot be removed or relocated after they have been released, their density can only increase over time at each point of the domain: the indeterminacy of the distribution of tracesterm·dom allows to satisfy this condition.

We demonstrate the validity of our framework in a variety of simulations in 1D and 2D. In 1D, we show that the swarmterm·key can achieve a stationary formation, and that we can make this stationary formation rigidly translate similar to a travelling waveterm·key in a non-dispersive medium. We validate our coordination strategy with even more complex scenarios, as we studied a concurrent modification of the formation during the translation of the wave, simulating a travelling waveterm·key in a dispersive medium, where different harmonics of the waves travel at a different speed, thus modifying the waveform over time. Analogous simulations at a discrete level demonstrated the validity and practical viability of our approach, even in very large-scale swarmsterm·key with thousands of units.

In 2D, we selected a historical example to show the potential of our approach. Specifically, we replicated the behaviour of Leonardo da Vinci’s lion-shaped automaton, which was designed to pay homage to the King of France through its complex movements. We showed that our approach allows the swarmterm·key to perform composite, time-modulated tasks, tracking an uneven trajectory, stopping at a precise location, modifying the formation in time, and resuming its motion along the trajectory. All these behaviours can be easily achieved through a stigmergic approach.

Our approach does not come without limitations. In particular, we acknowledge that further work is needed to support real-world implementation. The lack of an explicit model for the traceterm·dom deposition process benefits the generality of the theory, but it limits a prompt translation of the proposed approach to practice. Further, robustness of the approach to different forms of disturbances should be tested, and additional research is required for deployment in collective AMterm·dom. Further, we acknowledge that the proposed control strategy is not suitable for all collective tasks, especially for small teams in which coordination and localization within an absolute frame are more easily achieved. Nevertheless, we believe that the proposed mathematical framework offers a powerful and adaptable tool to design environmental modifications for swarmterm·key coordination, leaving designers the versatility to implement them in different ways.

Acknowledgements We are grateful to Mario di Bernardo for pointing us to the literature on stigmergyterm·dom and discussions related to continuificationterm·key, Chen Feng and Ludovic Righetti for helping us frame the collective constructionterm·key problem, Pietro De Lellis and Gian Carlo Maffettone for discussions related to the technical aspects of the work, and Petras Swissler for helping us position our work within the swarmterm·key robotics literature.

Contributor Information Alain BoldinipersName, Email: aboldini@nyit.edu.

Martina CivitellapersName, Email: mc8932@nyu.edu.

Maurizio PorfiripersName, Email: mporfiri@nyu.edu.

Ethics This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility All the data are included in the paper or in the code accessible at GitHub [ 72 ].

Electronic supplementary material is available online [ 73 ].

Declaration of AI use We have not used AI-assisted technologies in creating this article.

Authors’ contributions A.B.: conceptualization, formal analysis, investigation, methodology, software, visualization, writing—original draft, writing—review and editing; M.C.: formal analysis, methodology, software, writing—review and editing; M.P.: conceptualization, funding acquisition, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration We declare we have no competing interests.

Funding This work was supported by the National Science Foundation under grant no. CMMI 1932187 awarded to M.P.

References Associated Data Data Availability Statement All the data are included in the paper or in the code accessible at GitHub [ 72 ].

Electronic supplementary material is available online [ 73 ].