Case-Study: Boids Model - Modelling and Simulation Cookbook

The depict the concept of emergent behaviour via a simple agent-based model for flocking.

Topics:

Agent-Based Modelling

Objectives:

What does the term emergence mean in an agent-based model?
How does the runtime of agent-based model scale with the number of agents?

from collections import defaultdict

import numpy as np
from IPython.display import HTML
from matplotlib import animation
from matplotlib import pyplot as plt
from tqdm.notebook import tqdm

from modelling_and_simulation_cookbook import Clock

Background¶

In 1987, artificial-life and computer graphics expert Craig Reynolds published a groundbreaking model for the simulation and visualisation of flocking behaviour, the Boids model Reynolds (1987). His main motivation for model devleopment was that the complex motion of bird-flocks (e.g. starlings) or fish swarms (e.g. anchovies) is rarely seen in computer animation. That means, the overall purpose of the model is not to quantitatively recreate (or even predict) the behaviour of a flock, but to generate a simple realistic visualisation of one.

Conceptual Model¶

Agents¶

The model stands out for its simplicity and can be named as a prime example of a spatially continuous agent-based model. The $N$ agents in the model, also known as boids, represent the individual animals in the swarm and have two states: a position and a velocity. Let $i\in \{1,\dots,N\}$ refer to the $i$ -th agent in the model, then $\vec{p}_i(t)\in \mathbb{R}^{2}$ refers to the agent’s position at time $t$ and $\vec{v}_i(t)\in \mathbb{R}^{2}$ for the corresponding velocity. As a result, the whole swarm is given by a set of points and corresponding vecloity vectors - and can be visualised as such.

Time-Update General¶

The model is updated in discrete time steps whereas, velocities and positions of all boids are updated simultaneously, to make sure that no bias occurs due to iteration ordering.

Velocity Update - Interaction Radii¶

The most important part this update process is, how an individual boid $i$ , henceforth called the target boid, updates its velocity vector $\vec{v}_i(t)\rightarrow \vec{v}_i(t+1)$ . Key aspect for this is the defition of two radii, an observation radius $d_o\in \mathbb{R}^+$ and a collision radius $d_c\in \mathbb{R}^+$ , with $d_o>d_c$ . As seen in Figure 1, these two radii define two sets of other boids which are selected for interaction in the corresponding timestep for updating the velocity of the target boid (red):

I_d:=\{j\in\{1,\dots,i-1,i+1,\dots,N\}:\lVert\vec{p}_j(t)-\vec{p}_j(t)\rVert_2\leq d_o\},

(1)

I_c:=\{j\in\{1,\dots,i-1,i+1,\dots,N\}:\lVert\vec{p}_j(t)-\vec{p}_j(t)\rVert_2\leq d_c\}.

(2)

Due to $d_o>d_c$ clearly $I_c\subseteq I_o$ .

Definition of the boids selected for interaction in the given timestep via two radii d_o and d_c. — Figure 1:Definition of the boids selected for interaction in the given timestep via two radii $d_o$ and $d_c$ .

Velocity Update - Components¶

The update for the velocity is motivated by three general concepts: cohesion, alignment and separation.

Cohesion: Each boid aims to stay within the center of the swarm to stay safe. Therefore, the target boid is attracted to the centre of the boids within observation radius:
$\vec{w}_1:=\frac{1}{|I_o|}\sum_{j\in I_o}\vec{p}_j-\vec{p}_i.$
(3)
Alignment: Each boid aims to swim into the same direction as the swarm. Therefore, the target boid is attracted to the average velocity of the boids within observation radius:
$\vec{w}_2:=\frac{1}{|I_o|}\sum_{j\in I_o}\vec{v}_j.$
(4)
Separation: Finally, the boid does not want to run into its colleagues to avoid getting hurt. Therefore it is distracted from the centre of the boids within collision radius:
$\vec{w}_3:=-\frac{1}{|I_c|}\sum_{j\in I_o}\vec{p}_j-\vec{p}_i.$
(5)
These three components are visualised in Figure 2.

Although these three rules are, in principle, sufficient to generate meaningful swarm behaviour, it is often useful to add additional rules or vectors to make the update more elegant. We therefore add a retraction component

\vec{w}_0:= v_i(t)

(6)

to prevent the boid from changing its velocity too abruptly, and a memory component

\vec{w}_4:=-\vec{p}_i(t)\delta_{\lVert\vec{p}_i(t)\rVert_\infty>0.5}

(7)

to pull boids that have strayed too far from the origin back into the unit square $[-0.5,0.5]^2$ .

With these components, the velocity is updated via

\vec{v}_i^{raw}(t+1)=\vec{w}_0l_0+\vec{w}_1l_1+\vec{w}_2l_2+\vec{w}_3l_3+\vec{w}_4l_4,

(8)

whereas $l_0,\dots,l_4$ are cefficients between 0 and 1. Finally, we also consider a maximum velocity $v_{max}\in \mathbb{R}^+$ to truncate the velocities when necessary:

\vec{v}_i(t+1)=\vec{v}_i^{raw}(t+1)\frac{min(v_{max},\lVert \vec{v}_i^{raw}(t+1)\lVert_2)}{\lVert \vec{v}_i^{raw}(t+1)\lVert_2}.

(9)

Figure 2:Three main vector components for velocity-update in the boids model: cohesion, alignment, separation.

Note, that the extension of the model from two to three dimensions works without any further redefinition of the model.

Implemented Model¶

Model implementation can be made very straight forward and even independent of the space-dimension. The class below gives an example including routines for plotting.

class BoidsModel:
    def __init__(
        self,
        N: int,
        dim: int = 3,
        v_max: float | None = None,
        d_o: float | None = None,
        d_c: float | None = None,
        ls: list[float] | None = None,
        seed: int = 12345,
    ) -> None:
        """
        initialises a boids model
        :param N: number of boids
        :param dim: in how many dimensions the boids swim (usually 2 or 3)
        :param v_max: maximum velocity of the boids
        :param d_o: observation radius of the boids
        :param d_c: collision radius of the boids
        :param ls: vector of parameters for the velocity update. Components are:
          ls[0] ... maintain current direction,
          ls[1] ... direction of neighbor’s average position,
          ls[2] ... average direction of neighbor’s,
          ls[3] ... direction to avoid collision with other fish,
          ls[4] ... attraction to center
        :param seed: seed for the random number generator to make things reproducible
        """
        np.random.seed(seed)
        self.N = N
        self.dim = dim
        # select good default parametes
        if v_max is None:
            v_max = 0.03
        if d_o is None:
            d_o = 0.1
        if d_c is None:
            d_c = 0.04
        if ls is None:
            ls = [0.5, 0.1, 0.8, 0.8, 0.01]
        self.v_max = v_max  # maximal velocity of fish 3d
        self.d_o = d_o  # observation radius of fish in 3d
        self.d_c = d_c  # collision radius of fish 3d
        self.l0 = ls[0]  # ... maintain current direction
        self.l1 = ls[1]  # ... direction of neighbor’s average position
        self.l2 = ls[2]  # ... average direction of neighbor’s
        self.l3 = ls[3]  # ... direction to avoid collision with other fish
        self.l4 = ls[4]  # ... attraction to center

        self.v_max_sq = (
            self.v_max**2
        )  # compute square distance (so we do not need to compute roots)
        self.d_c_sq = (
            self.d_c**2
        )  # compute square distance (so we do not need to compute roots)
        self.d_o_sq = (
            self.d_o**2
        )  # compute square distance (so we do not need to compute roots)

        self.create_boids()  # initialise positions

    def create_boids(self) -> None:
        """
        Initialises the boids within the unit square/cube
        """

        self.pos = np.random.random([self.N, self.dim]) - 0.5
        self.vel = np.zeros([self.N, self.dim])

    def compute_rel_pos(
        self, own_pos: np.ndarray, nbs: np.ndarray | None = None
    ) -> tuple[np.ndarray, np.ndarray]:
        """
        Utility function to compute the relative position of a boid to an array
        of other boids
        :param own_pos: position of the target boid as 1 x dim vector
        :param nbs: optional array of indices to filter the original list
        :returns: relative positions to the boids (defined by the index array "nbs")
        as k x dim matrix and squared distances as k x 1 vector
        """
        nbs_pos = self.pos if nbs is None else self.pos[nbs, :]
        nbs_relpos = nbs_pos - own_pos[None, :]
        return nbs_relpos, (nbs_relpos * nbs_relpos).sum(axis=1)

    def compute_mean(
        self, array: np.ndarray, boolarray: np.ndarray | None = None
    ) -> np.ndarray:
        """
        Computes the mean of given array along the first axis.
        A boolarray can be given to filter the array.
        :param array: k x dim array or e.g. positions or velocities
        :param boolarray: k x 1 vector of boolean values to filter
        :returns: mean values as 1 x dim vector
        """
        if boolarray is None:
            return array.mean(axis=0)
        else:
            return array[boolarray, :].mean(axis=0)

    def update(self) -> None:
        """
        Performs an update step of the boids model
        """
        vel_new = np.zeros([self.N, self.dim])
        zero = np.zeros(self.dim)  # preallocate a zero vector to speed things up
        for i in range(self.N):
            own_pos = self.pos[i, :]
            own_vel = self.vel[i, :]
            w0 = own_vel
            relpos, distances = self.compute_rel_pos(own_pos)
            observed_bool = distances <= self.d_o_sq
            observed_bool[i] = False
            if observed_bool.any():
                w1 = self.compute_mean(relpos, observed_bool)
                w2 = self.compute_mean(self.vel, observed_bool)
                collison_bool = distances <= self.d_c_sq
                collison_bool[i] = False
                if collison_bool.any():
                    w3 = -self.compute_mean(relpos, collison_bool)
                else:
                    w3 = zero
            else:
                w1 = w2 = w3 = zero
            w4 = zero if not np.max(np.abs(own_pos)) > 0.5 else -own_pos
            veln = (
                self.l0 * w0 + self.l1 * w1 + self.l2 * w2 + self.l3 * w3 + self.l4 * w4
            )
            total_vel_sq = (veln * veln).sum()
            if total_vel_sq > self.v_max_sq:
                veln = veln * self.v_max / (total_vel_sq) ** 0.5
            vel_new[i, :] = veln
        self.vel = vel_new
        self.pos = self.pos + vel_new

    def plot_state_2d(self) -> None:
        """
        Creates a 2d quiver plot with matplotlib
        """
        cmp = plt.get_cmap("jet")  # slow = blue, fast = red
        total_vel_sq = (self.vel * self.vel).sum(axis=1)
        colors = [cmp(x / self.v_max_sq) for x in total_vel_sq]
        plt.scatter(
            self.pos[:, 0], self.pos[:, 1], s=1, c=colors
        )  # to also show those with 0-velocity
        plt.quiver(
            self.pos[:, 0],
            self.pos[:, 1],
            self.vel[:, 0],
            self.vel[:, 1],
            color=colors,
            scale=2.0,
            cmap="jet",
            clim=(0, self.v_max_sq),
        )

    def plot_state_3d(self) -> None:
        """
        Creates a 3d quiver plot with matplotlib
        """
        cmp = plt.get_cmap(
            "inferno"
        )  # things in the background are yellow, foreground almost black
        q = (
            self.pos[:, 0] - self.pos[:, 1] + self.pos[:, 2]
        )  # quantify what is in the foreground and what is in the background
        colors = [cmp((x + 1.5) / 3) for x in q]
        ax = plt.gca()
        ax.scatter(
            self.pos[:, 0], self.pos[:, 1], self.pos[:, 2], s=1, c=colors
        )  # to also show those with 0-velocity
        ax.quiver(
            self.pos[:, 0],
            self.pos[:, 1],
            self.pos[:, 2],
            self.vel[:, 0],
            self.vel[:, 1],
            self.vel[:, 2],
            colors=colors,
            normalize=False,
        )

Matplotlibs FuncAnimation provides a framework for animating the output. Unfortunately, it requires to compute all frames up-front. However, outside of an iPython notebook, also real-time animations are possible using different strategies.

def draw_frame(it: int, mdl: BoidsModel, dim: int) -> None:
    """
    Routine to be used in FuncAnimation to create a 2D or 3D video of
    the boids model
    :param it: iteration count
    :param mdl: boids model instance
    :param dim: dimension (2 or 3)
    """
    mdl.update()
    if dim == 2:
        plt.cla()
        mdl.plot_state_2d()
        plt.xlim([-0.7, 0.7])
        plt.ylim([-0.7, 0.7])
    elif dim == 3:
        ax = plt.gca()
        ax.cla()
        mdl.plot_state_3d()
        ax.set_xlim([-0.7, 0.7])
        ax.set_ylim([-0.7, 0.7])
        ax.set_zlim([-0.7, 0.7])


def animate_boids_model(
    mdl: BoidsModel, dim: int, frames: int, fps: float = 2.0
) -> HTML:
    """
    Creates a 2D or 3D animation of a boids model run over a certain
    number of update steps
    :param mdl: boids model instance
    :patam dim: dimension (2 or 3)
    :param frames: number of frames/update steps
    :param fps: frames per second for the animation
    :returns: html encoded video
    """
    plt.figure(figsize=(5, 5))
    if dim == 3:
        plt.gcf().add_subplot(projection="3d")
    plt.gcf().set_dpi(80)
    anim = animation.FuncAnimation(
        plt.gcf(),
        draw_frame,
        frames=tqdm(range(frames)),
        interval=1 / fps * 100,
        fargs=(mdl, dim),
    )
    video = HTML(anim.to_jshtml())
    plt.close()
    return video


# run and animate for 2D
mdl = BoidsModel(100, 2)
display(animate_boids_model(mdl, 2, 100))

The model clearly shows fascinating crowd behaviour just from applying the three (five) simple rules on the indiviual level. The macro behaviour is, up to some extent, completely unpredictable and therefore emergent. This makes the boids model the perfect case study to illustrate the corresponding key propery of agent-based models described by Bonabeau in 2002 Bonabeau (2002).

# run and animate for 3D
mdl = BoidsModel(100, 3, d_o=0.2)
display(animate_boids_model(mdl, 3, 100))

As seen, not much is going on in the 3D animation when using 100 boids. Clearly, we would like to increase the number of agents to replicate the size of real swarms. Note that, for example, a flock of starlings or a school of anchovies can easily consist of 50000 animals. We run a view computation time test to evaluate the performance.

def time_test(N: int):
    clock = Clock()
    mdl = BoidsModel(N, 2)
    clock.start()
    for _i in range(10):
        mdl.update()
    tm = clock.stop(f"Model with {N} agents")
    print(f"seconds/thousand_boids: {(tm / N * 1000).total_seconds()}")


time_test(256)
time_test(512)
time_test(1024)
time_test(2048)

Performance¶

The most important observation from the tests is not, that the computation time increases with $N$ - this had to be expected, because the update step consists mainly of a large loop over $N$ - but how it increases. Dividing the computation time by the number of boids $N$ still gives an increasing sequence, meaning that doubling the number of agents takes more than twice as long. In other words, the computation time increases nonlinearly with the number of agents which is one of the most problematic features of the boids model and one of the classic problems of agent-based models in general.

We spot the origin for this nonlinearity directly in the implementation: for every boid, the method compute_rel_pos computes the relative distances from the target boid to all other boids to determine which other boids are in the direct surrrounding. Although this computation is implemented well vectorised, it essentially translates to a loop over all other agents and therefore scales linearly with the number of boids. So, in the worst case, we experience an $\mathcal{O}(N^2)$ runtime, because for every boid ( $N$ ) we need to compute the distances to every other boid ( $N$ ).

Although a lot of performance can be gained by implementation tricks like vectorisation, this problem will always persist, unless one is willing to perform drastic changes. One idea is provided by tesselation as visualised in Figure 3. When dividing the whole area into equivalent zones of size $d_o\times d_0$ is becomes drastically easier to find boids in range $d_o$ . Suppose, the target boid is currently in zone $(i,j)$ , the number of boids which could be in observation zone are limited to the boids in the nine zones $\{(i+k,j+l),k\in \{-1,0,1\},l\in \{-1,0,1\}\}$ . That means, distance calculation only needs to be performed for the boids in these zones which significantly saves time, for the price of zone-assigment of the boids. An implementation is given below.

$Tesselation of the area in eqivalent zones of size d_o\times d_0. Instead of evaluating the distance to all other boids (left) to get the number of boids in range d_o, we only need to search in neighboured zones (right).$

Figure 3:Tesselation of the area in eqivalent zones of size $d_o\times d_0$ . Instead of evaluating the distance to all other boids (left) to get the number of boids in range $d_o$ , we only need to search in neighboured zones (right).

class BoidsModelZones(BoidsModel):
    def __init__(self, *args, **kwargs) -> None:
        """
        Performance improved implementation of the boids model using tesselation
        """
        super().__init__(*args, **kwargs)
        assert self.dim in {2, 3}, ValueError(
            "this implementation only works for 2D or 3D"
        )

    def create_boids(self) -> None:
        """
        Initialises the boids within the unit square/cube.
        Override the base method because it also creates a dictionary assigning zones to
        boids
        """
        super().create_boids()
        self.regions = np.floor(self.pos / self.d_o).astype(
            np.int32
        )  # rounds every position to the next multiple of d_o
        self.region_dict = defaultdict(
            set
        )  # creates a default dictionary with emtpy sets
        for i in range(self.N):
            reg = tuple(self.regions[i, :])  # convert to tuple to becom hashable
            self.region_dict[reg].add(i)

    def update_region_dict(self) -> None:
        """
        Updates the dictionary for assigment of zones to boids.
        This is faster than creating it new every time
        """
        new_regions = np.floor(self.pos / self.d_o).astype(np.int32)
        changed = np.any(
            new_regions != self.regions, axis=1
        )  # checks which regions have changed
        for i in np.nonzero(changed)[0]:
            old_reg = tuple(self.regions[i, :])
            new_reg = tuple(new_regions[i, :])
            self.region_dict[old_reg].remove(i)
            self.region_dict[new_reg].add(i)
        self.regions = new_regions

    def get_neighbors(self, i: int) -> np.ndarray:
        """
        Returns all indices to boids in the neighbouring 9 regions of the boid with
        index i
        """
        nbs = set()
        reg = self.regions[i, :]
        vc = (-1, 0, 1)
        if self.dim == 2:
            for j1 in vc:  # first coordinate
                for j2 in vc:  # second coordinate
                    ind = (reg[0] + j1, reg[1] + j2)
                    nbs.update(self.region_dict.get(ind, {}))
        else:
            for j1 in vc:  # first coordinate
                for j2 in vc:  # second coordinate
                    for j3 in vc:  # third coordinate
                        ind = (reg[0] + j1, reg[1] + j2, reg[2] + j3)
                        nbs.update(self.region_dict.get(ind, {}))
        nbs.remove(i)  # remove the target boid's index from the set
        return np.fromiter(
            nbs, dtype=np.int32
        )  # create a numpy array from the set for indexing

    def update(self):
        vel_new = np.zeros([self.N, self.dim])
        outs = (
            np.max(np.abs(self.pos), axis=1) > 0.5
        )  # vectorize comutation of who is outside of the observation area
        zero = np.zeros(self.dim)
        for i in range(self.N):
            own_pos = self.pos[i, :]
            own_vel = self.vel[i, :]
            nbs = self.get_neighbors(
                i
            )  # limit all computations to the boids in the neighboring zones
            w0 = own_vel
            if len(nbs) > 0:
                nbs_relpos, nbs_distances = self.compute_rel_pos(own_pos, nbs)
                observed_bool = nbs_distances <= self.d_o_sq
                # note that we do not need to manually set observed_bool[i]=False
                # because the target boid is not in "nbs" anyway
                if observed_bool.any():
                    w1 = self.compute_mean(nbs_relpos, observed_bool)
                    nbs_vel = self.vel[nbs, :]
                    w2 = self.compute_mean(nbs_vel, observed_bool)
                    collison_bool = nbs_distances <= self.d_c_sq
                    if collison_bool.any():
                        w3 = -self.compute_mean(nbs_relpos, collison_bool)
                    else:
                        w3 = zero
                else:
                    w1 = w2 = w3 = zero
            else:
                w1 = w2 = w3 = zero
            w4 = zero if not outs[i] else -own_pos
            vel_new[i, :] = (
                self.l0 * w0 + self.l1 * w1 + self.l2 * w2 + self.l3 * w3 + self.l4 * w4
            )

        total_vel_sq = (vel_new * vel_new).sum(
            axis=1
        )  # also vectorize computation of velocity truncating
        too_fast = total_vel_sq > self.v_max_sq
        facs = self.v_max / (total_vel_sq[too_fast]) ** 0.5
        vel_new[too_fast, :] = vel_new[too_fast, :] * facs[:, None]
        self.vel = vel_new
        self.pos = self.pos + vel_new
        self.update_region_dict()  # update zone assignment after movement

def time_test_zones(N: int):
    clock = Clock()
    mdl = BoidsModelZones(N, 2)
    clock.start()
    for _i in range(10):
        mdl.update()
    tm = clock.stop(f"Model with {N} agents")
    print(f"seconds/thousand_boids: {(tm / N * 1000).total_seconds()}")


time_test_zones(256)
time_test_zones(512)
time_test_zones(1024)
time_test_zones(2048)

The performance has clearly improved, but still does not scale linearly. This seems surprising, but is due to the increasing density of the boids. Increasing $N$ witout increasing the size of the grid (or decreasing the radii) at the same time causes more boids in each gridcell and consequently longer computation times. Nevertheless we are now able to perform comparably fast computations even with increased number of boids:

# run and animate for 3D
mdl = BoidsModelZones(500, 3, d_o=0.15, d_c=0.06)
display(animate_boids_model(mdl, 3, 100))

Parameter Studies¶

With the performance-improved implementation we can finally put our focus on what is actually going on in the model. In the following, we vary a few model parameters and have a look on what happens qualitatively with our swarm. As expected, the size of the observation and collision radii $d_o$ and $d_c$ are likely the most influencial parameters of the model next to the agent count $N$ .

d_o = 0.2
d_c = 0.1
print(f"Increased radii: {d_o}/{d_c}")
mdl = BoidsModelZones(100, 2, d_o=d_o, d_c=d_c)
display(animate_boids_model(mdl, 2, 150))
d_o = 0.04
d_c = 0.02
print(f"Decreased radii: {d_o}/{d_c}")
mdl = BoidsModelZones(100, 2, d_o=d_o, d_c=d_c)
display(animate_boids_model(mdl, 2, 150))

As seen, the radii must be chosen in a sensible range compared to the overall size of the observed area in order to produce realistic-looking swarm behaviour. If chosen too large, the swarm is stuck within the boundaries and has no space to move and evolve, if the values are too small, the resulting swarm(s) are tiny and tight and it will take very long until all boids are caught eventually.

We finally vary the weights $l_0,\dots,l_4$ .

ls = [0.5, 0.0, 0.8, 0.8, 0.01]
print(f"No cohesion: {ls}")
mdl = BoidsModelZones(100, 2, ls=ls)
display(animate_boids_model(mdl, 2, 100))
ls = [0.5, 0.1, 0.0, 0.8, 0.01]
print(f"No alignment: {ls}")
mdl = BoidsModelZones(100, 2, ls=ls)
display(animate_boids_model(mdl, 2, 100))
ls = [0.0, 0.1, 0.8, 0.8, 0.01]
print(f"No memory: {ls}")
mdl = BoidsModelZones(100, 2, ls=ls)
display(animate_boids_model(mdl, 2, 100))
ls = [0.5, 0.1, 0.8, 0.0, 0.01]
print(f"No separation: {ls}")
mdl = BoidsModelZones(100, 2, ls=ls)
display(animate_boids_model(mdl, 2, 100))
ls = [0.5, 0.8, 0.8, 0.8, 0.01]
print(f"High attraction to neighbors: {ls}")
mdl = BoidsModelZones(100, 2, ls=ls)
display(animate_boids_model(mdl, 2, 100))

First of all, we note that the three basic rules – cohesion, alignment and separation – do indeed represent minimum requirements, and that none of the three may be disregarded for the model. – Although the scenario without cohesion appears to produce sensible results, in this case the only thing holding the boids together is the repulsive force outside the observed area. This means that, were it not for this force, the swarm would slowly fall apart.

The results for the scenarios without alignment and without memory are quite similar. In both cases, there is no collective movement of the swarm in a specific direction. In the former case, each boid simply attempts to position itself within a safe zone between collision and observation, without aiming for a specific destination. In the latter case, any directional movement of the swarm is immediately halted as soon as a counteracting force arises.
The scenario without separation is very easy to interpret. As there are no dispersive forces at work, the positions within the swarm slowly converge due to cohesion. Consequently, a swarm in which all boids occupy the same position and have the same velocity represents a stable equilibrium state, which is reached if the simulation runs for a sufficiently long time. Such a trivial equilibrium state does not occur in any of the other scenarios.
The scenario in which the cohesion force between neighbours is increased shows very similar results to the scenario without separation. It can therefore be inferred that the cohesive force counteracts the repulsive force and that the interaction of the relevant parameters influences the observed distance between the swarm members.

References¶

Reynolds, C. W. (1987). Flocks, herds and schools: A distributed behavioral model. Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques, 25–34. 10.1145/37401.37406
Bonabeau, E. (2002). Agent-based modeling: Methods and techniques for simulating human systems. Proceedings of the National Academy of Sciences, 99(suppl_3), 7280–7287.