From ant colonies to termite systems to bee swarms, collective systems operate to solve complex problems efficiently. Somehow, it all happens without centralized control or direct communications between individual agents. Instead, the system relies on the concept of stigmergy, where agents modify the environment they work in, and those modifications alone become enough to guide the behavior of other agents. Through Agent Based Modeling (ABM) and cellular automata approaches, this project aims to implement a realistic model of stigmergic coordination, in which a population of agents aims to complete as many tasks as possible. Each agent’s movement is governed by two environmental factors: an attractive field derived from distances to the locations of tasks, encouraging agents to move towards tasks, and a repulsive field deposited from nearby agents left to diffuse and decay over time, discouraging grouping of agents. Each agent also possesses the ability of a short term memory, containing the most recent cells visited, which are forbidden to the agent, lessening repeated movements and increasing grid exploration. After exploring a fixed square grid, the model was extended to different shapes, supporting user chosen lattices for the agents’ regions, and was recreated in three dimensions, allowing real world situations to be modeled in irregular, constrained grids, rather than just idealized ones. An interactive interface allows the user to explore the governing parameters and explore the important relationships in the simulation and how local parameter choices can shape global coordination patterns. Moving beyond single run observations, several hundred simulations were done across varying parameter combinations, with results being graphed to explore how decay rate, memory length, and attraction/repulsion exponents can affect task completion time.
Introduction
Introduction
In 2017, researchers Danielli A. Lima and Gina M.B. Oliveira published their novel paper, A cellular automata ant memory model of foraging in a swarm of robots, detailing a theoretical approach to using Agent Based Modeling (ABM) to simulate stigmergy, a process by which a population of agents works together to complete tasks without any form of central communication. The approach the two scholars outlined relied on the idea of “inverted pheromones”, where agents deposit pheromones that repel other agents, in an attempt to explore the search area more efficiently. Given the Wolfram Language’s specialty in cellular automata, this bagged the question: could this approach to modeling the phenomenon of stigmergy be carried out in Wolfram? This research project aims to answer that question by combining the theory of “inverted pheromones” with the classical approach of “attractive tasks”, modeling how the agents explore a search grid filled with tasks, and establishing relationships through data collection and repeated simulation. Beyond simply confirming that such a model could be constructed in the Wolfram Language, this project extends the approach to arbitrary geometries and characterizes, across thousands of simulation runs, how the model’s governing parameters shape the emergent coordination behavior of the swarm.
Random Movement
Random Movement
First, before implementing attracted to robots and agents repulsed from each other, I must model the simplest kind of agent movement: randomness. At each tick, each agent will make an unweighted random choice between moving up, left, right, and down, confined to a 20x20 grid, and if it happens to move onto a task, that task will be counted as “completed” and deleted.
To start off, I define a moveAgents function that can take a coordinate and randomly pick one of the four cells in its Von Neumann neighborhood to move to.
Defines a function to make a random choice for moving up, left, down, or right given a coordinate:
In[]:=
moveAgentRandomly[{x_,y_},gridSize_]:=Module[{listChoices,validChoices},listChoices={x,y}+#&/@({{0,-1},{0,1},{1,0},{-1,0}});validChoices=DeleteCases[listChoices,{a_,b_}/;(a<1||a>gridSize||b<1||b>gridSize)];RandomChoice[validChoices]]
Now, because I will represent the agents as a list of coordinates, I must to have some way to apply this to a list. Of course, I cannot have two agents moving onto the same coordinate, so I will also add a safeguard against that.
Defines a function that randomly moves a list of agents and blocks two agents from moving to the same square:
In[]:=
moveMultipleAgentsRandomly[agents_,gridSize_]:=Module[{result},result=moveAgentRandomly[#,gridSize]&/@agents;If[DuplicateFreeQ[result],result,moveMultipleAgentsRandomly[agents,gridSize]]]
Building on this, I can package this all into a function that can take in a lists of agents and tasks, move the agents, and delete any tasks that an agent moved on to.
Sets up a function to continuously update the agents and tasks when called. Every time the agents are moved, any tasks sharing a coordinate with an agent are deleted, and if the tasks list becomes empty, it is reset to a list of random coordinates:
In[]:=
stepStateRandomly[{agents_,tasks_},gridSize_,numTasks_,validCoordinates_]:=Module[{newAgents,newTasks},newAgents=moveMultipleAgentsRandomly[agents,gridSize];newTasks=Select[tasks,!MemberQ[newAgents,#]&];If[newTasks=={},newTasks=RandomChoice[validCoordinates,numTasks]];{newAgents,newTasks}]
Before visualizing the output, I will write a function that can take a list of agents and a list of tasks, and display them effectively.
In[]:=
frameGraphicsRandom[{agents_,tasks_},gridSize_]:=Graphics[Join[{Red},Disk[#,0.5]&/@agents,{Green},Disk[#,0.5]&/@tasks],PlotRange->{{0,gridSize},{0,gridSize}}];
Subsequently, I can visualize the result of the random movement. Note that randomness is not very effective at reaching the tasks, as expected.
Displays the output of continuously randomly moving the agents. Note the use of NestList to generate a list of the first 100 steps of moving the agents, and the slider which controls which step is currently being displayed:
In[]:=
Manipulate[DynamicModule[{gridSize=20,validCoordinates,initialAgents,initialTasks,stateFrames},validCoordinates=Flatten[Table[{i,j},{i,1,gridSize},{j,1,gridSize}],1];initialAgents=RandomChoice[validCoordinates,numAgents];initialTasks=RandomChoice[validCoordinates,numTasks];stateFrames=NestList[stepStateRandomly[#,gridSize,numTasks,validCoordinates]&,{initialAgents,initialTasks},100];Dynamic[frameGraphicsRandom[stateFrames[[currentState]],gridSize]]],{{numAgents,5,"Number of Agents"},Range[10],ControlType->Slider,Appearance->"Labeled"},{{numTasks,5,"Number of Tasks"},Range[10],ControlType->Slider,Appearance->"Labeled"},{{currentState,1,"Current Step"},Range[100],ControlType->Slider,Appearance->"Labeled"},SaveDefinitions->True,TrackedSymbols:>{numAgents,numTasks,currentState}]
Tabu Memory
Tabu Memory
Having seen how inefficient simple random movement is, let’s try something to optimize the patterns of the agents a little bit more: Tabu Memory. This algorithm works by letting each agent possess a memory containing its most recent steps. If an agent is next to a coordinate in its memory, it is not allowed to move onto that cell, even if the cell contains a task. I can implement this by changing the structure of how I store the positions of agents: instead of just a list of coordinates, I use a list of associations, where each agent corresponds to an association. In each association will be a coordinate, for location, and a list of past coordinates, for the memory.
First, I will create a new function to move agents that will use queues and not allow agents to move to spots in their memory. When an agent has to move, it still makes a random choice between up, left, right, and down, but, after adding Tabu memory, it is not allowed to move into a cell in its memory. After an agent moves into a cell, that cell is added to its memory, and if the size of the memory is equal to the maximum size allowed by the code (tabuMemoryLength), the first element of the memory is removed.
Defines another agent movement function that takes into account that each agent is now an association. Note how originalListChoices is a list of the four adjacent cells to the agent, validChoices is a list of cells that are valid within the 20x20 grid, and memoryValidChoices is a list of choices that are both within the grid and not in the memory:
In[]:=
moveAgentTabu[agent_,gridSize_,tabuMemoryLength_]:=Module[{originalListChoices,validChoices,memoryValidChoices,newPos,newQueue},originalListChoices=agent["pos"]+#&/@({{0,-1},{0,1},{1,0},{-1,0}});validChoices=DeleteCases[originalListChoices,{a_,b_}/;(a<1||a>gridSize||b<1||b>gridSize)];memoryValidChoices=Select[validChoices,!MemberQ[Normal[agent["Tabu"]],#]&];If[memoryValidChoices=={},newPos=RandomChoice[validChoices],newPos=RandomChoice[memoryValidChoices]];newQueue=agent["Tabu"];If[Length[Normal[newQueue]]==tabuMemoryLength,newQueue["Pop"],""];newQueue["Push",newPos];<|"pos"->newPos,"Tabu"->newQueue|>]
I will also define a new function that moves multiple agents and does not allow two agents into the same cell, accounting for the fact that each agent is represented by an association.
Defines another function to move multiple agents, again disallowing two agents moving to the same coordinate on a tick. The variable agents is still a list of associations, but note how I access all the coordinates using a pure anonymous function on the result:
In[]:=
moveMultipleAgentsTabu[agents_,gridSize_,tabuMemoryLength_]:=(Module[{result},result=moveAgentTabu[#,gridSize,tabuMemoryLength]&/@agents;If[DuplicateFreeQ[#["pos"]&/@result],result,moveAgentTabu[#,gridSize,tabuMemoryLength]&/@agents]])
I define another function to update the positions of agents and tasks that is compatible with the tabu memory. Similar to when I used completely random movement, at each tick, only two things happen: each agent is given a new position, and any task sharing a position with an agent is deleted.
Declares a function to move the agents whenever called, keeping in mind that each agent is an association:
In[]:=
stepStateTabu[{agents_,tasks_},numTasks_,gridSize_,tabuMemoryLength_,validCoordinates_]:=Module[{newAgents,newTasks,listAgentCoords},newAgents=moveMultipleAgentsTabu[agents,gridSize,tabuMemoryLength];listAgentCoords=#["pos"]&/@agents;newTasks=Select[tasks,!MemberQ[listAgentCoords,#]&];If[newTasks=={},newTasks=RandomChoice[validCoordinates,numTasks],""];{newAgents,newTasks}]
Because the format of how the agents are stored has changed, I must also update how the visualization is constructed.
Similar to stepStateTabu, in this function declaration to create a simple visualization of the agents and tasks, the list of the coordinates of the agents, the list of coordinates is found by applying a pure anonymous function to agents:
In[]:=
frameGraphicsTabu[{agents_,tasks_},gridSize_]:=(Module[{listAgentCoords},listAgentCoords=#["pos"]&/@agents;Graphics[Join[{Red},Disk[#,0.5]&/@listAgentCoords,{Green},Disk[#,0.5]&/@tasks],PlotRange->{{0,gridSize},{0,gridSize}}]])
Finally, I visualize the results again. I can see that, in contrast to the example above, the agents do not double back on themselves at all, leading to a more efficient grid search.
Visualizes our tabu memory implementation. Again, I use NestList to get a list of all frames, and control which frame is currently being used with the currentStep variable:
In[]:=
Manipulate[DynamicModule[{gridSize=20,validCoordinates,initialAgents,initialTasks,stateFrames},validCoordinates=Flatten[Table[{i,j},{i,1,gridSize},{j,1,gridSize}],1];initialAgents=Table[<|"pos"->RandomChoice[validCoordinates],"Tabu"->CreateDataStructure["Queue"]|>,numAgents];initialTasks=RandomChoice[validCoordinates,numTasks];stateFrames=NestList[stepStateTabu[#,numTasks,gridSize,tabuMemoryLength,validCoordinates]&,{initialAgents,initialTasks},100];Dynamic[frameGraphicsTabu[stateFrames[[currentState]],gridSize]]],{{numAgents,5,"Number of Agents"},Range[10],ControlType->Slider,Appearance->"Labeled"},{{numTasks,5,"Number of Tasks"},Range[10],ControlType->Slider,Appearance->"Labeled"},{{currentState,1,"Current Step"},Range[100],ControlType->Slider,Appearance->"Labeled"},{{tabuMemoryLength,4,"Length of Tabu Memory"},Range[10],ControlType->Slider,Appearance->"Labeled"},SaveDefinitions->True,TrackedSymbols:>{numAgents,numTasks,tabuMemoryLength}]
Repulsive Agents
Repulsive Agents
However, despite my addition of memory, the agents are still extremely inefficient, and for good reason: they have no reason to explore the graph in a timely manner. Let’s try adding something else, in addition to memory: repulsive agents, where each agent repels the others by dropping repulsive pheromones. Using this technique, each agent will be motivated to explore different sections of the search area from its counterparts, thus making the grid search more efficient, as tasks will be found quicker.
To start implementing the concept of repulsive agents, first, I will implement a function to move the agents, taking into account the repulsive values. I can accomplish this by having a matrix that is the same size as a search grid, called repulsiveVals, where each cell’s value is its corresponding repulsive value. Then, when an agent has to move, it observes the repulsive values of the four adjacent cells, and instead of making a completely random choice to decide where to go, it will make a weighted random choice using the repulsive values as weights. Consequently, the agents will begin avoiding each other, making grid search much more efficient.
I define a new function to move the agents, taking the repulsive weights into account. Observe the addition of the variables repulsiveVals, repulsiveExponent, and weights, and note that when there are valid cells to move on, newPos is determined by a weighted random choice. See how repulsiveExponent controls how intense the attractiveness of a cell is - the higher repulsiveExponent is, the more likely an agent is to to a cell with a higher repulsive value. Also, note that (1/1+repulsiveVals) is used, because the lesser the actual value of the matrix is, the more attractive it is:
In[]:=
moveAgentRepulsive[agent_,gridSize_,tabuMemoryLength_,repulsiveVals_,repulsiveExponent_]:=Module[{originalListChoices,validChoices,memoryValidChoices,newPos,newQueue,weights},originalListChoices=agent["pos"]+#&/@({{0,-1},{0,1},{1,0},{-1,0}});validChoices=DeleteCases[originalListChoices,{a_,b_}/;(a<1||a>gridSize||b<1||b>gridSize)];memoryValidChoices=Select[validChoices,!MemberQ[Normal[agent["Tabu"]],#]&];weights=(1/(1+repulsiveVals[[#[[1]],#[[2]]]]))^(repulsiveExponent)&/@memoryValidChoices;If[memoryValidChoices=={},newPos=RandomChoice[validChoices],newPos=RandomChoice[weights->memoryValidChoices]];newQueue=agent["Tabu"];If[Length[Normal[newQueue]]==tabuMemoryLength,newQueue["Pop"],""];newQueue["Push",newPos];<|"pos"->newPos,"Tabu"->newQueue|>]
Similarly, I must update the function to move multiple agents, taking the repulsive values into account.
Builds a function to move multiple agents simultaneously. The repulsiveVals matrix has become a new function argument, as it is needed to move the agents:
In[]:=
moveMultipleAgentsRepulsive[agents_,gridSize_,tabuMemoryLength_,repulsiveVals_,repulsiveExponent_]:=(Module[{result},result=moveAgentRepulsive[#,gridSize,tabuMemoryLength,repulsiveVals,repulsiveExponent]&/@agents;If[DuplicateFreeQ[#["pos"]&/@result],result,moveAgentRepulsive[#,gridSize,tabuMemoryLength,repulsiveVals,repulsiveExponent]&/@agents]])
Despite all this, the code to frame the agents and tasks does not change.
In[]:=
frameGraphicsRepulsive[{agents_,tasks_},gridSize_]:=(Module[{listAgentCoords},listAgentCoords=#["pos"]&/@agents;Graphics[Join[{Red},Disk[#,0.5]&/@listAgentCoords,{Green},Disk[#,0.5]&/@tasks],PlotRange->{{0,gridSize},{0,gridSize}}]])
Because of all these changes, I must update the recursive function to step between agent states. I also add the idea of pheromone decay here, which entails multiplying the matrix of repulsive values by a constant less than one each time to ensure that as the agent moves away, the spot is no longer repulsive.
Constructs a function to step through the agents’ movements, implementing the agents’ repulsiveness towards each other:
Now that the agents have the functionality to react to the cells around them being repulsive, I must actually implement the idea that the cells are repulsive. I will do this by having, every tick, each agent add to the repulsive values of the cells in its Moore neighborhood (the 3x3 matrix centered at that cell). Through the pheromone decay, the repulsiveness of the cells will slowly fade to zero, making it look like each agent leaves a trail of repulsiveness.
To put this idea into action, first, I must be able to get the Moore neighborhood of a coordinate, and validate that its cells are inside our 20x20 grid.
Takes in a coordinate and returns its Moore neighborhood, the 3x3 matrix centered at that coordinate. The function automatically deletes any cells outside the 20x20 grid:
Here is my function to update the repulsive values matrix. Given a coordinate, addMatrixRepulsive will update the values of the repulsive matrix, by adding to the repulsive values of its Moore neighborhood.
Takes in the location of a single agent and updates the repulsiveVals matrix accordingly. The critical part is using the offsets such as #[[1]]-(coordinate[[1]]-2) to get the right term to add:
Analogous to the move agents function, I must create a function that can apply this to all agents, which is what updateRepulsiveVals does. After applying updateRepulsiveVals to a list of coordinates (the positions of the matrix) and the old repulsive matrix, a new repulsive matrix will be returned containing updated values
Takes in a list of all agents’ positions, and performs the update to the repulsiveVals grid for all of them:
I also define the function that simulates decay of pheromones, multiplying each value in the repulsive matrix by a constant less than one.
Defines a function to take a matrix and multiply all of its values by a number less than one, effectively simulating decay:
I can again visualize the results of the repulsive values simulation. Note how the agents seem to avoid each other more and more. Included is a heat map detailing the relative probabilities that an agent moves to a square, where the darker a square is, the less likely it is for an agent to move there. As expected, each agent seems to have a “trail” of dark squares in the heat map, whose length can be varied by adjusting the pheromone decay rate.
Visualizes the agents and their tasks, along with a heat map detailing probabilities:
Attractive Tasks
Attractive Tasks
There is some improvement in the efficiency of how the agents explore the grid after implementing both Tabu Memory and repulsive features, but for the tasks to be completed as fast as possible, I need to complete the final piece of the puzzle: making the tasks attractive to the agents.
This can be done by using a system similar to the repulsive values: giving each coordinate an “attractiveness value” based on how far the Manhattan Distance is to the nearest task. Just like with the repulsive values, I use the attractive values as weights, enabling the agents to make weighted random choices, becoming more likely to move closer to the tasks and completing them more efficiently.
First, I define the function to move the agents using the attractive values as weights. It is analogous to the functions for the repulsive weighted agents, but with an attractive values matrix.
Encodes a function to move agents by having them make weighted random choices based off of the attractive values of the surrounding cells. Observe how, similar to the repulsive case, I employ (1/(1+attractiveVals)) because the lower the value of the attractive values matrix is, the more “attractive” that cell is, as it is closer to a task. Also, notice how attractiveExponent controls the intensity of the attractiveness - the higher it is, the more likely an agent is to move to a cell with a better attractive value than to move to a cell with a worse attractive value:
Just as before, I must also define a function to apply this to more than one agent.
Creates a function building off the last one, which moves multiple agents simultaneously:
My function to plot the agents and disks does not need to be changed for this application, as the structure of how the agents and tasks are stored remains the same.
Defines a function to visualize the agents and tasks analogous to those that came before it:
Now, I address the problem of actually updating the attractive values of each cell. I can do this by treating the grid as a graph, where each cell is a node connected by an edge to the cells above, to the left, to the right, and below it. Then, I can efficiently use the plethora of Wolfram’s built in graph searching functions to find the Manhattan distance from each cell to its nearest task.
To begin, I encode the functionality to take in a coordinate and the graphical representation of the grid, and find the Manhattan distance to the nearest task.
Defines a function that can take in a coordinate, the graph of the coordinate grid, and the tasks, and can output the Manhattan distance to the nearest task, or that coordinate’s attractive value:
Then, I create a function that applies this to all cells of the coordinate grid I work within, which becomes my method of updating the attractive vals.
Creates a function that finds the attractive value of every single coordinate. Note how I normalize each vector, because most values will be much bigger than 1:
Having completed those functions, I can write a function to step through moving the agents, analogous to the one for the repulsive agents.
Declares a function to step through moving the agents, taking attractive values into account. At each step, three things happen: the agents are updated, making weighted choices based off the attractive values, tasks that match coordinates of agents are deleted, and if the tasks have changed, then the attractive values are updated. Observe how the attractive values only need to be changed if the tasks are changed:
Now, I visualize the results of the attractive agents logic with a heat map.
Again, this code is relatively simple given all the defined functions. It initializes the variables, creates a list of all frames through a NestList, and enables a slider to move through the frames and see the heat map progress:
The Combination
The Combination
Now that I have implemented both repulsive agents and attractive tasks separately, it is time to take the two methods, and combine them to fully optimize the way the agents explore the grid. As all the necessary functions are already almost entirely defined, it is only a matter of updating the functions that move the agents and step through the grid to fully optimize the agents. To make matters simpler, I will employ a third matrix, called totalVals, that houses each cell’s total value: the combination of its repulsiveness and attractiveness, and only use that in deciding where an agent moves.
Creates a new function to move agents that accounts for attractive and repulsive factors. Every time an agent moves, it looks at the total value of the cells adjacent to it, and makes a weighted random choice using those values as weight:
As before, I must define a function to apply the new function for moving agents to all agents, with a safeguard against two agents moving on to the same cell.
Declares a function to move multiple agents simultaneously, with each agent using totalVals as its weights:
Lastly, I must define the function to step through each agent movement, which will generate frames as it goes. This is the most complex state stepping function yet, and each tick, it does multiple things. After the agents and tasks are updated, the function uses applyPheromoneDecay and updateRepulsiveVals to update the repulsive values matrix, and if it is confirmed that the tasks have changed, the attractive values matrix is updated. Then, the new total value is calculated for each cell, using my (1/(1+x)) logic, and the repulsive and attractive exponents are used to tune the intensity of the values. The list of new agents, tasks, attractive values, repulsive values, and total values is returned, and in the later calculation, this will be repeated many times.
As always, at the end, I put the Manipulate and Dynamic Module, allowing me to play with the parameters and explore the relationships between them interactively.
Custom Shapes
Custom Shapes
Now that I have implemented both repulsive agents and attractive tasks on an idealized square in a two dimensional grid, the natural direction to go from here is to extend this to other shapes on the coordinate plane. The ability to model stigmergy on the irregular shapes of nature instead of an unrealistic square will unlock many new capabilities, making it possible to study coordination in settings that actually resemble the environments these systems are meant to operate in. To incorporate this, I have the user specify a list of lattice points before running the function, and use the minimal convex region encapsulating all those points as the confining region for the agents and tasks.
To implement this new feature, I need to remake all of the previous functions, as more complexity has been added into the model.
First, given the user’s list of coordinates, these functions find the smallest convex region containing all of them, as the user could have entered the coordinates in a mistaken order, leading to a difficult, possibly self intersecting, shape to work with.
Defines two functions that, given the list of coordinates the user inputs, can output the minimal convex region containing all the points and the boundary of that region:
Having done that, I have to redefine all of the old functions to account for this added complexity, because the user may have fed in coordinates that are offset from the origin, leading to out of bounds error when trying to feed those coordinates in alone as indices of a matrix.
Establishes a new function to move a singular agent that takes into account that the coordinates the user fed in may be offset from the origin:
Also, in the new visualization, I include the user’s polygon at the very back layer to emphasize that the agents and tasks are confined to that region.
Constructs a new function to visualize the agents and tasks at each step. Here, the convex polygon the user has specified is placed at the bottom layer of the graphics to emphasize that the agents and tasks are confined to that region:
I must also construct a new function to enable moving multiple agents simultaneously, given this added complexity. Note how the bottom left corner of the rectangle has been included as a function input to account for user coordinate offset from the origin.
Builds a new function to move multiple agents at once that is compatible with the moveAgentCustom function:
I redefine the old mooreNeighborhood function here for convenience, which can take a coordinate and output the 3x3 matrix centered at that coordinate.
Redefines the old function to get the Moore neighborhood of a coordinate:
Similarly to the function to move the agents, I must redefine my functions to update the repulsive and attractive values to account for coordinates that may be offset from the origin.
Updates the functions that changed the repulsive values to account for offset:
Here, the logic for updating the attractive values stays largely the same: treat the coordinate grid as a graph, where each coordinate is a node connected by an edge to the coordinates above, to the left, to the right, and below it. Because I use Table instead of directly accessing the matrix’s elements through indexing, there is no need to account for offset here.
Creates new functions to update the attractive values of the cells, where the logic and inputs stay largely unchanged compared to other functions:
I also create a new function to update the total values matrix. Most of the updates here come from accounting for offset from the origin.
Declares a function to update the total values matrix, accounting for offset from the origin in the user’s input:
My final function is the classic stepState function, which has significantly more arguments than previous editions, stemming from the added complexity of user coordinate input. At each tick, the familiar process happens: agents and tasks are updated, the attractive values are updated if the tasks have changed, the repulsive values are updated with the new values of the agents, and the total values are updated based off the attractive and repulsive values of each cell.
Creates a function that takes in all the input of the user combined with the parameters necessary, and updates all the values one tick at a time:
Finally, at last, I can visualize the results of the custom shape feature. Using an approximation of a regular hexagon, I can see that the agents and tasks are confined to the bounds of the shape, and successfully only move within it.
Visualizes the output of the custom coordinates feature using an approximation of a hexagon of side length 10 :
Three Dimensions
Three Dimensions
Extending the framework from two dimensions to three is a natural next step. The underlying logic of agent movement, task seeking, and stigmergic signaling carries over largely unchanged. However, what does change is the environment itself, which now approximates real-world spatial structure rather than a simplified plane. To make this transition, I re-implement each of the original functions, adapting them to operate over a third axis of movement.
Defines the moveAgent function with a three dimensional adaption. Note the addition of the withinCube function and the cubeSize variable, as it becomes tedious to continuously manually check if a point exists within a cube:
I must also make a new function to frame the agent and tasks. Observe how they are now spheres instead of disks:
The aforementioned withinCube function, to easily check if a point is within the given cube by simply making sure each of the three x, y, and z coordinates are valid:
Here, I update the familiar function to find the three dimensional Moore neighborhood of a point. This is achieved by constructing the three by three by three cuboid lattice of points centered at that point, and then deleting any points outside the valid region:
The following three functions are the adaptations of the repulsive agents logic, extended to three dimensions.
Defines a function to update the repulsive values of the surrounding locations of an agent in three dimensions:
Encodes a function responsible for decaying the repulsive values of a matrix. The code remains exactly the same as the previous functions to decay the repulsive values, but I redefine it here to avoid confusion:
As before, updateRepulsiveVals takes in a list of agents, and updates the repulsive values of each cell in order to push agents away from each other:
Finally, I encode the functions to update the attractiveness of each task as well as the total values for each location.
Given a coordinate, the movement space as a graph, and a list of tasks, findDistanceToTasks3D will find the manhattan distance from that coordinate to the nearest task:
Creates a function responsible for updating the attractive value of each coordinate when called. Note how it is necessary to normalize the values because most of them will be far bigger than 1:
The final piece of the puzzle, this function will update the combined attractive and repulsive value of each coordinate. Observe how the intensity of the values can be tuned by adjusting attractiveMultiplier and repulsiveMultiplier:
Having completed all of the minor functions necessary, the last thing to do is to combine them all into a recursive function that will, as the name suggests, step through each of the states.
Implements a function to step through moving the agents and the tasks.
Lastly, I visualize the results, using the function NestList as before to generate a list of frames and displaying each progressive frame with a slider.
Parameter Optimization
Parameter Optimization
With the simulation fully implemented across both two and three dimensions, the model is now ready to turn from being built to being evaluated. Rather than treating the parameter choices used so far as fixed, I run a series of single parameter sweeps, varying one of the region topology (square, hexagon, cube), attractive exponent, repulsive exponent, Tabu memory length, or pheromone decay rate at a time while keeping all other parameters constant. The goal is twofold: to find parameter regimes that improve performance, and to look for qualitative patterns in the resulting plots that might reveal something about the underlying stigmergic dynamics.
First, I construct functions that measure the performance of a certain set of values for the parameters. Here, I will use how long the five agents take to complete five tasks as our metric.
This function, stepStateTiming, is the recursive function I will use to step through moving agents one at a time and measure how long it takes for all five tasks to be deleted . It has the same syntax as the earlier function stepStateCustom, but with one crucial change . In the previous function, for display purposes, whenever there were zero tasks remaining, the tasks would be reset randomly . However, as the goal here is to realize when there are no tasks left and stop, that line must be removed .
This function, calcRunTimeOnce2D, allows me to feed in any list of points, and it will return the result of how long it took 5 agents to complete 5 tasks in one simulation .
Defines a baseline function that, given any list of points, simulates the agents moving and returns an experimental value for how long it took for all tasks to be completed. Observe how NestWhile is used to continuously step the agents until the list of tasks is empty, at which point the amount of time spent to reach all five tasks is returned.
Now, I can define functions for a square and a hexagon that, for a given set of values for the attractive exponent, repulsive exponent, pheromone decay rate, and tabu memory length, will return the result of how long it took to reach all five tasks in some number of simulations.
Defines a function that, specifically for the 20x20 square grid, can return the result of any number of simulations with the parameter values specified in the input.
Defines a similar function that, specifically for a hexagon approximating a regular one of side length 10, can return the result of any number of simulations with the parameter values specified in the input.
I use the same idea with a 3D model, creating functions to run simulations and return how long it took all tasks to be completed.
Again, I must use an altered version of our previous function to step through the three dimensional agents’ movements, as the tasks should not automatically be reset upon all being reached:
Creates a baseline function for running one simulation inside a 3D shape and returning how long it took for all tasks to disappear:
Analogous to its 2D counterparts, this function will take in the four parameters, and return the results of running any number of simulations within a cube and seeing how long it took for all tasks to disappear .
Defines a high level function to return the result of running any number of simulations using the parameters specified in the input:
For each of the four parameters of interest, I can keep the other 3 parameters the same and vary the fourth across a plausible range of its optimal values, and see what trends the graph produces. However, I did choose to vary the attractive exponent and the repulsive exponent together, given their symbiotic relationship.
Here is an example of how I can produce a graph I can visually analyze, using the method to vary the pheromone decay rate of the agents on a two dimensional square grid. First, I can set the values of the attractive exponent, repulsive exponent, and tabu memory length to constants, and vary the pheromone decay rate. For each value of the pheromone decay rate, I will run the simulation 50 times, and take the mean of all the found runtimes.
Using that, I calculate a list of coordinates for this dataset.
However, only calculating a list of coordinates for one dataset is not very reliable. To rectify this, we run this simulation for four other sets of values for the attractive exponent, the repulsive exponent, the decay rate, and the tabu memory length.
Calculates data points by running the simulation for four other parameter sets. For example, dataset3 sets the attractive exponent to 40, the repulsive exponent to 40 and the Tabu Memory length to 8, and varies the pheromone decay rate:
Turns the data I have generated into coordinates I can plot on a graph:
Using the same method for all of our parameters of interest, I generated the following 9 plots, each with 4-5 different datasets to help observe trends.
Attractive Exponent vs. Repulsive Exponent vs. Mean Time Spent on a 7x7x7 Cube:
Attractive Exponent vs. Repulsive Exponent vs. Mean Time Spent on a Regular Hexagon Side Length 10:
Attractive Exponent vs. Repulsive Exponent vs. Mean Time Spent on a 20x20 Square Grid:
Pheromone Decay Rate vs. Mean Time Spent on a 20x20 Grid:
Pheromone Decay Rate vs. Mean Time Spent on a 7x7x7 Cube:
Pheromone Decay Rate vs. Mean Time Spent on a Regular Hexagon Side Length 10:
Tabu Memory Length vs. Mean Time Spent on a 7x7x7 Cube:
Tabu Memory Length vs. Mean Time Spent on a 20x20 Grid:
Tabu Memory Length vs. Mean Time Spent on a Regular Hexagon Side Length 10:
Through analysis of these graphs, there are a few patterns I noticed that can be explored in further studies. First, throughout all the graphs with attractive exponents, repulsive exponents, and the mean time spent, past a certain threshold for the repulsive exponent, the timeSpent seems oddly linear in terms of the repulsive exponent. Second, it seems to be a general trend that as the length of the Tabu Memory increases, the mean run time decreases, a pattern worth investigating. Third, oddly, it does not seem that the mean run time follows a consistent run time as a function of the pheromone decay rate, which seems out of place to me, as how quickly a trail decays directly controls how much agents can control a grid.
Through more detailed investigation and rigorous statistical testing, I am sure that many more patterns can be discovered, paving a path for the field of studying stigmergy.
Through more detailed investigation and rigorous statistical testing, I am sure that many more patterns can be discovered, paving a path for the field of studying stigmergy.
Conclusion
Conclusion
In this project, I implemented a thorough Agent Based Model of stigmergic coordination, where a population of agents aims to complete as many tasks as possible without any centralized communication or direct control. Using cellular automata, I implemented many different features of the agents on a 2D square to explore what led to the most efficient grid search, such as memory, repulsiveness of agents, and attractiveness of tasks. Then, after fully exploring a simply 20x20 grid, I moved to exploring less idealized shapes, implementing custom lattices for the agents’ bounding regions, and even recreating the model in a three dimensional 7x7x7 cube, allowing me to model the nuances of everyday life. Through interactive interfaces, data collection and graphical analysis, I was able to identify some patterns that could be investigated in the future, which is itself the more interesting takeaway: a handful of simple local rules produce global patterns complex enough that they aren’t obvious even to someone who wrote the rules, the hallmark of cellular automata.
Limitations
Limitations
The primary limitation of this research project came with the data collections and analysis. Despite dedicating days of computing to the procedure, I could only run a few hundred simulations, and thus, there is not nearly enough evidence to back up my claims. Furthermore, all patterns that I see are simply based upon observation from the naked eye. I did not have the time or the knowledge to perform tests for statistical significance, and thus, cannot formally state anything as a certainty.
Future Directions
Future Directions
Building on top of the work that I have done here, there are many extensions that could be carried out to explore stigmergy within Agent Based Modeling. One weak spot in my model was the lack of obstacles to the agents, as in the real world, agents will not be allowed to move wherever they please. To help solve this problem, one could implement something like a stationary or moving predator: if an agent moves onto a spot containing a predator, the agent is deleted. Another weakness of my program was the inefficiency of the code, with multiple DynamicModules wrapped inside of Manipulates, something that could easily be improved upon while preserving the quality of the program.
Furthermore, a promising direction one could extend this project is towards doing more research into the relationship between the parameters and environment, and the performance of the model. Maybe researchers could look more into how changing between differently shaped environments would affect how quickly tasks could be completed. Perhaps they could create a general case formula, writing an expression for the expected mean time for any number of tasks to be completed in terms of the number of agents, the tabu memory length, the attractive and repulsive exponents, and the pheromone decay rate. Or, researchers could investigate other metrics for how successful the grid search was, such has how many times agents crossed each others paths, or how many turns an agent made. The world of stigmergy is just starting to open up with the computational tools available today, and I would love to see it taken further than what I have done.
Furthermore, a promising direction one could extend this project is towards doing more research into the relationship between the parameters and environment, and the performance of the model. Maybe researchers could look more into how changing between differently shaped environments would affect how quickly tasks could be completed. Perhaps they could create a general case formula, writing an expression for the expected mean time for any number of tasks to be completed in terms of the number of agents, the tabu memory length, the attractive and repulsive exponents, and the pheromone decay rate. Or, researchers could investigate other metrics for how successful the grid search was, such has how many times agents crossed each others paths, or how many turns an agent made. The world of stigmergy is just starting to open up with the computational tools available today, and I would love to see it taken further than what I have done.
Acknowledgements
Acknowledgements
I would like to thank my mentor Felipe Amorim for his continued support, expertise, and unending patience as I worked on this project. I would also like to thank Stephen Wolfram, Rory Foulger, Eryn Gillam, Megan Davis, Cyrus Taylor, Daniel Xu, Eric Archerman, Arnay Garhyan, and all the wonderful staff at the Wolfram Summer Research Program for making this opportunity possible for me. Huge thanks goes out to my advocate group: Suvid Bordia, Lukas Tsui, and Theo Archer, and to Jordana Wilkes for helping me out with some timely words of wisdom.
References
References
◼
Lima, D. A., & Oliveira, G. M. B. (2017). A cellular automata ant memory model of foraging in a swarm of robots [Review of A cellular automata ant memory model of foraging in a swarm of robots]. Applied Mathematical Modelling, 47, 551–572. https://doi.org/10.1016/j.apm.2017.03.021
◼
Coordination of robot swarms in crowded and confined. (n.d.). https://crablab.gatech.edu/pages/publications/pdf/AINA-DISSERTATION-2024.pdf
◼
Wilkes, J. (2026). [WSRP25] Modeling slime mold intelligence through cellular automata - Online Technical Discussion Groups—Wolfram Community [Review of [WSRP25] Modeling slime mold intelligence through cellular automata - Online Technical Discussion Groups—Wolfram Community]. Wolfram.Com. https://community.wolfram.com/groups/-/m/t/3500639
◼
Wei, H. (2026). [WSRP25] Modelling pedestrian dynamics and crowd flow with 2D cellular automata - Online Technical Discussion Groups—Wolfram Community [Review of [WSRP25] Modelling pedestrian dynamics and crowd flow with 2D cellular automata - Online Technical Discussion Groups—Wolfram Community]. Wolfram.Com. https://community.wolfram.com/groups/-/m/t/3502839
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
A cellular automata model of stigmergic agent swarms
by Aarav Gupta
Wolfram Community, STAFF PICKS, July 9, 2026
https://community.wolfram.com/groups/-/m/t/3752559
by Aarav Gupta
Wolfram Community, STAFF PICKS, July 9, 2026
https://community.wolfram.com/groups/-/m/t/3752559