Swarm robotics uses large numbers of relatively simple robots that interact locally to produce global behaviors. In this lesson we focus on three canonical objectives:

• Aggregation: agents converge to a tight cluster.
• Dispersion: agents spread out to maintain separation.
• Coverage: agents distribute to optimally monitor an area.

We model a swarm of \( N \) identical point robots in \( \mathbb{R}^d \) with first-order kinematics \( \dot{x}_i(t) = u_i(t) \) where \( x_i(t)\in\mathbb{R}^d \) is position and \( u_i(t)\in\mathbb{R}^d \) is the control input of robot \( i \). Robots sense neighbors through a fixed or distance-limited interaction graph.

Let \( G = (V,E) \) be an undirected graph with \( |V| = N \), adjacency matrix \( A=[a_{ij}] \) with \( a_{ij}=1 \) if robots \( i \) and \( j \) interact, and degree matrix \( D=\mathrm{diag}(d_i) \) with \( d_i = \sum_j a_{ij} \). The (combinatorial) Laplacian is

\[ L = D - A. \]

Stacking all positions as \( z = [x_1^\top,\dots,x_N^\top]^\top\in\mathbb{R}^{Nd} \), interaction laws of the form \( u_i = f_i(\{x_j - x_i: j\in\mathcal{N}_i\}) \) (with neighbor set \( \mathcal{N}_i \)) define global dynamics \( \dot{z} = F(z) \).

The design problem is: choose \( f_i \) so that the emergent steady-state configuration realizes aggregation, dispersion, or coverage, while respecting locality and symmetry constraints.

flowchart TD
  A["Specify swarm task: 'aggregate', 'disperse', or 'cover'"] --> B["Choose interaction graph G and sensing radius R"]
  B --> C["Design local law u_i = f_i(neighbor positions)"]
  C --> D["Analyze stability / convergence using Lyapunov or gradient methods"]
  D --> E["Implement discrete-time update x_i(k+1) = x_i(k) + h u_i(k)"]
  E --> F["Simulate and tune gains, neighborhood rules"]
        

A standard aggregation law is consensus (agreement on a common position). For each agent,

\[ \dot{x}_i = -k \sum_{j\in\mathcal{N}_i} (x_i - x_j), \quad k > 0. \]

In vector form, using the Kronecker product \( \otimes \) and identity \( I_d \),

\[ \dot{z} = -k (L \otimes I_d)\, z. \]

2.1 Invariance of the swarm average

Define the average position \( \bar{x}(t) = \frac{1}{N}\sum_{i=1}^N x_i(t) \). Differentiating and using the consensus law:

\[ \dot{\bar{x}}(t) = \frac{1}{N}\sum_{i=1}^N \dot{x}_i = -\frac{k}{N}\sum_{i=1}^N \sum_{j\in\mathcal{N}_i} (x_i - x_j). \]

Since the graph is undirected, each term \( x_i - x_j \) appears once as \( (i,j) \) and once as \( (j,i) \), and they cancel pairwise:

\[ \sum_{i=1}^N \sum_{j\in\mathcal{N}_i} (x_i - x_j) = 0 \quad \Rightarrow \quad \dot{\bar{x}}(t) = 0. \]

Thus the swarm average remains constant in time: \( \bar{x}(t) = \bar{x}(0) \). This is an important invariant: under consensus aggregation the swarm collapses to its initial center of mass.

2.2 Lyapunov analysis of aggregation

Let \( z\in\mathbb{R}^{Nd} \) be the stacked state and define the quadratic Lyapunov function

\[ V(z) = \frac{1}{2} z^\top (L\otimes I_d)\,z = \frac{1}{4}\sum_{i,j=1}^N a_{ij}\,\|x_i - x_j\|^2. \]

Because \( L \) is symmetric positive semidefinite for an undirected graph, we have \( V(z)\ge 0 \), and \( V(z)=0 \) if and only if \( x_1 = \dots = x_N \), i.e. perfect aggregation.

Using \( \dot{z} = -k(L\otimes I_d)z \), the derivative is

\[ \dot{V}(z) = z^\top (L\otimes I_d)\dot{z} = -k\, z^\top (L^2\otimes I_d)\, z \le 0, \]

since \( L^2 \) is also symmetric positive semidefinite. Therefore \( V(z) \) is non-increasing along trajectories.

The largest invariant set where \( \dot{V} = 0 \) is

\[ \mathcal{M} = \{z : (L\otimes I_d)\,z = 0\}. \]

For a connected graph, \( \ker L = \mathrm{span}\{ \mathbf{1} \} \), where \( \mathbf{1} \in \mathbb{R}^N \) is the all-ones vector. Hence,

\[ \mathcal{M} = \{\mathbf{1}\otimes x^\star : x^\star \in \mathbb{R}^d\}, \]

which is the consensus subspace. By LaSalle's invariance principle, all trajectories converge to \( \mathcal{M} \), and by average invariance the consensus value is \( x^\star = \bar{x}(0) \). Thus, aggregation to the initial average position is globally asymptotically stable.

Aggregation collapses robots into a point; dispersion instead enforces non-zero inter-agent distances, often around a desired spacing \( d_{\mathrm{des}} > 0 \). A common approach is to define a pairwise potential depending only on distance \( r_{ij} = \|x_i - x_j\| \):

\[ \Phi(r) = \frac{1}{4}(r^2 - d_{\mathrm{des}}^2)^2. \]

The total swarm potential is

\[ W(x_1,\dots,x_N) = \sum_{(i,j)\in E} \Phi(\|x_i - x_j\|), \]

and the control law is a gradient descent on this potential:

\[ \dot{x}_i = - \nabla_{x_i} W = -\sum_{j\in\mathcal{N}_i} \Phi'(\|x_i - x_j\|) \frac{x_i - x_j}{\|x_i - x_j\|}. \]

For our choice of \( \Phi \),

\[ \Phi'(r) = (r^2 - d_{\mathrm{des}}^2)\, r, \]

giving

\[ \dot{x}_i = -\sum_{j\in\mathcal{N}_i} (r_{ij}^2 - d_{\mathrm{des}}^2)\, (x_i - x_j). \]

If \( r_{ij} > d_{\mathrm{des}} \) then the term is attractive (robots move closer); if \( r_{ij} < d_{\mathrm{des}} \) then it is repulsive (robots move apart). Equilibria correspond to critical points of \( W \); configurations where all neighboring distances equal \( d_{\mathrm{des}} \) are stationary points.

3.1 Lyapunov function and stability

Consider the Lyapunov candidate equal to the potential, \( V = W \). Its time derivative along trajectories is

\[ \dot{V} = \sum_{i=1}^N \nabla_{x_i} W^\top \dot{x}_i = -\sum_{i=1}^N \|\nabla_{x_i} W\|^2 \le 0, \]

which implies \( W \) is non-increasing. Under mild regularity and connectivity conditions, the swarm converges to the set of configurations where all gradients vanish, i.e. shapes with all neighbor distances at critical values of \( \Phi \). By choosing \( \Phi \) with a single minimum at \( r = d_{\mathrm{des}} \) on the relevant range of distances, one enforces dispersion around that spacing.

In practice, one often adds saturation, damping, or obstacle potentials to ensure bounded velocities and collision avoidance with the environment.

Coverage is a distributed sensing problem: robots must position themselves so that a domain \( Q\subset\mathbb{R}^2 \) is well monitored, possibly according to a density \( \phi(q)\ge 0 \) describing importance or event likelihood. Let robot positions be \( P = (p_1,\dots,p_N) \).

4.1 Voronoi partition

The Voronoi cell of robot \( i \) is

\[ V_i(P) = \{ q\in Q : \|q - p_i\| \le \|q - p_j\|\ \forall j\neq i \}. \]

A standard quadratic coverage cost is

\[ H(P) = \sum_{i=1}^N \int_{V_i(P)} \|q - p_i\|^2 \phi(q)\,dq. \]

Define the mass and centroid of cell \( V_i \):

\[ m_i(P) = \int_{V_i(P)} \phi(q)\,dq, \quad c_i(P) = \frac{1}{m_i(P)} \int_{V_i(P)} q\,\phi(q)\,dq. \]

Using calculus of variations one can show

\[ \nabla_{p_i} H(P) = 2\, m_i(P)\, (p_i - c_i(P)). \]

4.2 Gradient descent and Lloyd's algorithm

A continuous-time coverage controller is

\[ \dot{p}_i = -k_c\, \nabla_{p_i} H(P) = -2k_c\, m_i(P)\, (p_i - c_i(P)). \]

The time derivative of \( H \) along trajectories is

\[ \dot{H}(P) = \sum_{i=1}^N \nabla_{p_i} H^\top \dot{p}_i = -4k_c \sum_{i=1}^N m_i(P)\, \|p_i - c_i(P)\|^2 \le 0. \]

Hence \( H \) decreases monotonically, and equilibria satisfy \( p_i = c_i(P) \) for all \( i \), i.e. the robots form a centroidal Voronoi tessellation.

A discrete-time implementation (Lloyd's algorithm) alternates:

1. Compute Voronoi cells \( V_i(P^k) \).
2. Compute centroids \( c_i(P^k) \).
3. Update \( P^{k+1} \leftarrow (c_1(P^k),\dots,c_N(P^k)) \).

This defines a descent method for \( H \), converging to centroidal configurations under mild conditions.

For numerical simulation or digital control, we discretize the swarm dynamics with step size \( h > 0 \):

\[ x_i(k+1) = x_i(k) + h\, u_i(k). \]

For consensus aggregation,

\[ x_i(k+1) = x_i(k) - h k_a \sum_{j\in\mathcal{N}_i} (x_i(k) - x_j(k)). \]

In vector form,

\[ z(k+1) = \bigl(I_{Nd} - h k_a (L\otimes I_d)\bigr)\, z(k). \]

A standard result: if \( 0 < h k_a < \frac{2}{\lambda_{\max}(L)} \), where \( \lambda_{\max}(L) \) is the largest eigenvalue of \( L \), then the discrete-time consensus dynamics are stable and converge to the consensus subspace.

A generic simulation loop (for any local law) can be summarized as:

flowchart TD
  S["Initialize positions x_i(0) and parameters"] --> L1["For each time step k"]
  L1 --> L2["For each robot i: find neighbors within radius R"]
  L2 --> L3["Compute control u_i(k) for aggregation / dispersion / coverage"]
  L3 --> L4["Update x_i(k+1) = x_i(k) + h u_i(k)"]
  L4 --> L5["Optionally add noise, saturation, and visualize"]
  L5 --> L1
        

We implement a simple 2D swarm with three selectable behaviors: aggregation, dispersion, and coverage (via approximate Lloyd iterations using random samples of the domain).

import numpy as np

class Swarm2D:
    def __init__(self, N=30, dt=0.05, behavior="aggregate",
                 d_des=0.2, domain=((0.0, 1.0), (0.0, 1.0))):
        self.N = N
        self.dt = dt
        self.behavior = behavior
        self.d_des = d_des
        self.domain = domain
        # positions: shape (N, 2)
        self.x = np.random.rand(N, 2)
        self.x[:, 0] = self.x[:, 0] * (domain[0][1] - domain[0][0]) + domain[0][0]
        self.x[:, 1] = self.x[:, 1] * (domain[1][1] - domain[1][0]) + domain[1][0]
        # parameters
        self.k_agg = 1.0
        self.k_disp = 0.5
        self.k_cov = 1.0
        self.R = 0.4  # sensing radius

    def neighbors(self, i):
        diffs = self.x - self.x[i]
        dists = np.linalg.norm(diffs, axis=1)
        idx = np.where((dists > 0.0) & (dists < self.R))[0]
        return idx, diffs[idx], dists[idx]

    def step_aggregate(self):
        u = np.zeros_like(self.x)
        for i in range(self.N):
            idx, diffs, _ = self.neighbors(i)
            if idx.size > 0:
                u[i] = -self.k_agg * np.sum(diffs, axis=0)
        self.x += self.dt * u

    def step_dispersion(self):
        u = np.zeros_like(self.x)
        for i in range(self.N):
            idx, diffs, dists = self.neighbors(i)
            if idx.size == 0:
                continue
            # Pairwise potential derivative (r^2 - d_des^2) r
            r2 = dists ** 2
            phi_prime = (r2 - self.d_des ** 2) * dists
            # unit vectors
            dirs = diffs / dists[:, None]
            # sum over neighbors
            force = -np.sum(phi_prime[:, None] * dirs, axis=0)
            u[i] = self.k_disp * force
        self.x += self.dt * u

    def step_coverage(self, M_samples=500):
        # Approximate Lloyd step: sample domain, assign to nearest robot,
        # compute sample-based centroids, move robots toward them.
        xs = np.random.rand(M_samples, 2)
        xs[:, 0] = xs[:, 0] * (self.domain[0][1] - self.domain[0][0]) + self.domain[0][0]
        xs[:, 1] = xs[:, 1] * (self.domain[1][1] - self.domain[1][0]) + self.domain[1][0]
        # density phi(q) can be non-uniform; here uniform
        # assign samples to closest robot
        dists = np.linalg.norm(xs[:, None, :] - self.x[None, :, :], axis=2)  # M x N
        closest = np.argmin(dists, axis=1)
        centroids = np.zeros_like(self.x)
        counts = np.zeros(self.N, dtype=int)
        for m in range(M_samples):
            i = closest[m]
            centroids[i] += xs[m]
            counts[i] += 1
        for i in range(self.N):
            if counts[i] > 0:
                c_i = centroids[i] / counts[i]
                # gradient step toward centroid
                self.x[i] += self.dt * self.k_cov * (c_i - self.x[i])

    def step(self):
        if self.behavior == "aggregate":
            self.step_aggregate()
        elif self.behavior == "disperse":
            self.step_dispersion()
        elif self.behavior == "cover":
            self.step_coverage()
        # Optional: keep robots inside domain (simple clipping)
        self.x[:, 0] = np.clip(self.x[:, 0], self.domain[0][0], self.domain[0][1])
        self.x[:, 1] = np.clip(self.x[:, 1], self.domain[1][0], self.domain[1][1])

if __name__ == "__main__":
    import matplotlib.pyplot as plt

    swarm = Swarm2D(N=40, behavior="aggregate")
    for k in range(200):
        swarm.step()
        if k % 20 == 0:
            plt.clf()
            plt.scatter(swarm.x[:, 0], swarm.x[:, 1])
            plt.xlim(swarm.domain[0])
            plt.ylim(swarm.domain[1])
            plt.title(f"Step {k}")
            plt.pause(0.01)
    plt.show()
      

Switching behavior between "aggregate", "disperse", and "cover" demonstrates the three swarm behaviors introduced in this lesson.

Below is a compact C++11 implementation of aggregation and dispersion for 2D point robots. For real robots, this logic would be embedded in each robot's local controller or in a simulation environment (e.g., ROS node).

#include <vector>
#include <cmath>
#include <random>

struct Vec2 {
    double x{0.0}, y{0.0};
    Vec2() = default;
    Vec2(double x_, double y_) : x(x_), y(y_) {}
    Vec2 operator+(const Vec2& o) const { return Vec2(x + o.x, y + o.y); }
    Vec2 operator-(const Vec2& o) const { return Vec2(x - o.x, y - o.y); }
    Vec2 operator*(double s) const { return Vec2(x * s, y * s); }
    Vec2& operator+=(const Vec2& o) { x += o.x; y += o.y; return *this; }
};

double norm(const Vec2& v) {
    return std::sqrt(v.x * v.x + v.y * v.y);
}

class Swarm2D {
public:
    enum Behavior { AGGREGATE, DISPERSE };

    Swarm2D(std::size_t N, double dt, Behavior b)
        : N_(N), dt_(dt), behavior_(b), x_(N) {
        std::mt19937 gen(42);
        std::uniform_real_distribution<double> uni(0.0, 1.0);
        for (auto& p : x_) {
            p.x = uni(gen);
            p.y = uni(gen);
        }
    }

    void step() {
        if (behavior_ == AGGREGATE) stepAggregate();
        else stepDispersion();
    }

    const std::vector<Vec2>& positions() const { return x_; }

private:
    std::size_t N_;
    double dt_;
    Behavior behavior_;
    std::vector<Vec2> x_;
    double R_ = 0.4;
    double k_agg_ = 1.0;
    double k_disp_ = 0.5;
    double d_des_ = 0.2;

    void neighbors(std::size_t i,
                   std::vector<std::size_t>& idx,
                   std::vector<Vec2>& diffs,
                   std::vector<double>& dists) const {
        idx.clear();
        diffs.clear();
        dists.clear();
        for (std::size_t j = 0; j < N_; ++j) {
            if (j == i) continue;
            Vec2 d = x_[j] - x_[i];
            double r = norm(d);
            if (r > 0.0 && r < R_) {
                idx.push_back(j);
                diffs.push_back(d);
                dists.push_back(r);
            }
        }
    }

    void stepAggregate() {
        std::vector<Vec2> u(N_);
        std::vector<std::size_t> idx;
        std::vector<Vec2> diffs;
        std::vector<double> dists;
        for (std::size_t i = 0; i < N_; ++i) {
            neighbors(i, idx, diffs, dists);
            Vec2 ui(0.0, 0.0);
            for (const auto& d : diffs) ui += d;
            u[i] = ui * (-k_agg_);
        }
        for (std::size_t i = 0; i < N_; ++i) {
            x_[i] += u[i] * dt_;
        }
    }

    void stepDispersion() {
        std::vector<Vec2> u(N_);
        std::vector<std::size_t> idx;
        std::vector<Vec2> diffs;
        std::vector<double> dists;
        for (std::size_t i = 0; i < N_; ++i) {
            neighbors(i, idx, diffs, dists);
            Vec2 ui(0.0, 0.0);
            for (std::size_t k = 0; k < idx.size(); ++k) {
                const Vec2& d = diffs[k];
                double r = dists[k];
                double r2 = r * r;
                double phi_prime = (r2 - d_des_ * d_des_) * r;
                Vec2 dir(d.x / r, d.y / r);
                ui += dir * (-phi_prime);
            }
            u[i] = ui * k_disp_;
        }
        for (std::size_t i = 0; i < N_; ++i) {
            x_[i] += u[i] * dt_;
        }
    }
};
      

Rendering and interaction with real-time robot middleware (ROS, etc.) can wrap this class to update positions and publish control commands.

The following Java code implements aggregation and dispersion with a simple object-oriented structure. It can be extended with GUI drawing or network communication to robots.

public class Vec2 {
    public double x, y;
    public Vec2(double x, double y) { this.x = x; this.y = y; }
    public Vec2 plus(Vec2 o) { return new Vec2(x + o.x, y + o.y); }
    public Vec2 minus(Vec2 o) { return new Vec2(x - o.x, y - o.y); }
    public Vec2 times(double s) { return new Vec2(x * s, y * s); }
    public double norm() { return Math.sqrt(x * x + y * y); }
}

public class Swarm2D {
    public enum Behavior { AGGREGATE, DISPERSE }
    private int N;
    private double dt;
    private Behavior behavior;
    private Vec2[] x;
    private double R = 0.4;
    private double kAgg = 1.0;
    private double kDisp = 0.5;
    private double dDes = 0.2;

    public Swarm2D(int N, double dt, Behavior b) {
        this.N = N; this.dt = dt; this.behavior = b;
        this.x = new Vec2[N];
        java.util.Random rand = new java.util.Random(42);
        for (int i = 0; i < N; ++i) {
            x[i] = new Vec2(rand.nextDouble(), rand.nextDouble());
        }
    }

    private void neighbors(int i,
                           java.util.List<Integer> idx,
                           java.util.List<Vec2> diffs,
                           java.util.List<Double> dists) {
        idx.clear(); diffs.clear(); dists.clear();
        for (int j = 0; j < N; ++j) {
            if (j == i) continue;
            Vec2 d = x[j].minus(x[i]);
            double r = d.norm();
            if (r > 0.0 && r < R) {
                idx.add(j);
                diffs.add(d);
                dists.add(r);
            }
        }
    }

    public void step() {
        switch (behavior) {
            case AGGREGATE: stepAggregate(); break;
            case DISPERSE: stepDispersion(); break;
        }
    }

    private void stepAggregate() {
        Vec2[] u = new Vec2[N];
        java.util.List<Integer> idx = new java.util.ArrayList<>();
        java.util.List<Vec2> diffs = new java.util.ArrayList<>();
        java.util.List<Double> dists = new java.util.ArrayList<>();
        for (int i = 0; i < N; ++i) {
            neighbors(i, idx, diffs, dists);
            Vec2 ui = new Vec2(0.0, 0.0);
            for (Vec2 d : diffs) ui = ui.plus(d);
            u[i] = ui.times(-kAgg);
        }
        for (int i = 0; i < N; ++i) {
            x[i] = x[i].plus(u[i].times(dt));
        }
    }

    private void stepDispersion() {
        Vec2[] u = new Vec2[N];
        java.util.List<Integer> idx = new java.util.ArrayList<>();
        java.util.List<Vec2> diffs = new java.util.ArrayList<>();
        java.util.List<Double> dists = new java.util.ArrayList<>();
        for (int i = 0; i < N; ++i) {
            neighbors(i, idx, diffs, dists);
            Vec2 ui = new Vec2(0.0, 0.0);
            for (int k = 0; k < idx.size(); ++k) {
                Vec2 d = diffs.get(k);
                double r = dists.get(k);
                double r2 = r * r;
                double phiPrime = (r2 - dDes * dDes) * r;
                Vec2 dir = new Vec2(d.x / r, d.y / r);
                ui = ui.plus(dir.times(-phiPrime));
            }
            u[i] = ui.times(kDisp);
        }
        for (int i = 0; i < N; ++i) {
            x[i] = x[i].plus(u[i].times(dt));
        }
    }

    public Vec2[] getPositions() { return x; }
}
      

A simple Swing or JavaFX application can periodically call step() and render getPositions().

MATLAB lends itself naturally to matrix-based swarm simulations. The following script implements aggregation and dispersion; a Simulink model can replicate the update rule with vector blocks and a discrete-time integrator.

function swarm_sim(behavior)
% behavior: 'aggregate' or 'disperse'

if nargin < 1, behavior = 'aggregate'; end

N = 40;
dt = 0.05;
T = 20.0;
K = round(T / dt);
R = 0.4;
k_agg = 1.0;
k_disp = 0.5;
d_des = 0.2;

x = rand(N, 2);  % positions

for k = 1:K
    u = zeros(N, 2);
    for i = 1:N
        diffs = x - x(i, :);
        dists = sqrt(sum(diffs.^2, 2));
        mask = (dists > 0) & (dists < R);
        nbr_diffs = diffs(mask, :);
        nbr_dists = dists(mask);
        if isempty(nbr_dists), continue; end
        switch behavior
            case 'aggregate'
                u(i, :) = -k_agg * sum(nbr_diffs, 1);
            case 'disperse'
                r2 = nbr_dists .^ 2;
                phi_prime = (r2 - d_des^2) .* nbr_dists;
                dirs = nbr_diffs ./ nbr_dists;
                force = -sum(phi_prime .* dirs, 1);
                u(i, :) = k_disp * force;
        end
    end
    x = x + dt * u;
    if mod(k, 10) == 0
        clf;
        scatter(x(:, 1), x(:, 2), 'filled');
        axis([0 1 0 1]); axis equal;
        title(sprintf('k = %d', k));
        drawnow;
    end
end
end
      

Simulink implementation sketch: represent each robot's dynamics as \( x_i(k+1) = x_i(k) + h\,u_i(k) \) using a Discrete-Time Integrator block. A MATLAB Function block can take the vector of all positions and output control inputs according to the chosen aggregation or dispersion law. For coverage, one can approximate centroids using numerical integration or sampling inside another MATLAB Function block.

Mathematica's concise list operations are convenient for prototyping swarm dynamics. Below we implement aggregation and dispersion as pure functions and simulate with NDSolve-style explicit updates.

ClearAll[aggregateStep, disperseStep];

neighbors[xs_, i_, R_] := Module[{xi, diffs, dists, mask},
  xi = xs[[i]];
  diffs = xs - xi;
  dists = Norm /@ diffs;
  mask = MapIndexed[#1 > 0 && #1 < R && First[#2] != i &, dists];
  Pick[Transpose[{Range[Length[xs]], diffs, dists}], mask]
];

aggregateStep[xs_, dt_, kAgg_, R_] := Module[{u},
  u = Table[
    Module[{nbrs, sumDiff = {0., 0.}},
      nbrs = neighbors[xs, i, R];
      If[nbrs =!= {},
        sumDiff = Total[nbrs[[All, 2]]];
      ];
      -kAgg * sumDiff
    ],
    {i, Length[xs]}
  ];
  xs + dt * u
];

disperseStep[xs_, dt_, kDisp_, R_, dDes_] := Module[{u},
  u = Table[
    Module[{nbrs, force = {0., 0.}},
      nbrs = neighbors[xs, i, R];
      If[nbrs =!= {},
        force = Total[
          Module[{d = #2, r = #3, r2, phiPrime, dir},
            r2 = r^2;
            phiPrime = (r2 - dDes^2) * r;
            dir = d / r;
            -phiPrime * dir
          ] & @@@ nbrs
        ];
      ];
      kDisp * force
    ],
    {i, Length[xs]}
  ];
  xs + dt * u
];

(* Simulation *)
Nrobots = 30;
dt = 0.05;
steps = 400;
R = 0.4;
kAgg = 1.0;
kDisp = 0.5;
dDes = 0.2;

(* Initial positions *)
xs = RandomReal[{0, 1}, {Nrobots, 2}];

behavior = "aggregate"; (* or "disperse" *)

traj = NestList[
  If[behavior === "aggregate",
    aggregateStep[#, dt, kAgg, R] &,
    disperseStep[#, dt, kDisp, R, dDes] &
  ],
  xs, steps
];

ListAnimate[
  ListPlot[#, PlotRange -> { {0, 1}, {0, 1} }, AspectRatio -> 1,
    PlotStyle -> Directive[PointSize[0.02]]] & /@ traj
]
      

Extending this to coverage control requires computing approximate cell centroids, e.g. by sampling points in the domain and using Nearest to assign samples to robots.

Problem 1 (Average Invariance in Consensus): Consider the continuous-time consensus law \( \dot{x}_i = -k \sum_{j\in\mathcal{N}_i} (x_i - x_j) \) on a connected undirected graph. Prove that the average position \( \bar{x}(t) = \frac{1}{N}\sum_{i=1}^N x_i(t) \) is constant in time.

Solution:

Compute the time derivative:

\[ \dot{\bar{x}}(t) = \frac{1}{N}\sum_{i=1}^N \dot{x}_i = -\frac{k}{N}\sum_{i=1}^N\sum_{j\in\mathcal{N}_i}(x_i - x_j). \]

For an undirected graph, \( (i,j)\in E \) if and only if \( (j,i)\in E \). Rewriting the double sum over edges,

\[ \sum_{i=1}^N\sum_{j\in\mathcal{N}_i}(x_i - x_j) = \sum_{(i,j)\in E} \bigl[(x_i - x_j) + (x_j - x_i)\bigr] = 0. \]

Therefore \( \dot{\bar{x}}(t)=0 \), implying \( \bar{x}(t) = \bar{x}(0) \) for all \( t \).

Problem 2 (Lyapunov Function for Aggregation): For the same consensus system, define \( V(z) = \frac{1}{4}\sum_{i,j} a_{ij}\|x_i - x_j\|^2 \). Show that \( \dot{V} \le 0 \) and that \( V(z)=0 \) if and only if \( x_1 = \dots = x_N \).

Solution:

Using \( L = D - A \), \( V(z) = \frac{1}{2} z^\top (L\otimes I_d) z \). Then

\[ \dot{V} = z^\top (L\otimes I_d)\dot{z} = z^\top (L\otimes I_d)\bigl(-k(L\otimes I_d)z\bigr) = -k\, z^\top (L^2\otimes I_d) z \le 0, \]

since \( L^2 \) is symmetric positive semidefinite. The Laplacian of a connected graph has a one-dimensional nullspace spanned by \( \mathbf{1} \), so \( V(z)=0 \) if and only if \( (L\otimes I_d)z = 0 \), i.e. all \( x_i \) are equal.

Problem 3 (Coverage Gradient): Let \( Q\subset\mathbb{R}^2 \), \( \phi(q)\ge 0 \), and coverage cost \( H(P) = \sum_{i=1}^N \int_{V_i(P)} \|q - p_i\|^2 \phi(q)\,dq \) . Derive \( \nabla_{p_i} H(P) \) under the assumption that Voronoi cell boundaries move negligibly with \( p_i \).

Solution:

Neglecting boundary motion, the derivative is

\[ \frac{\partial H}{\partial p_i} = \int_{V_i(P)} \frac{\partial}{\partial p_i} \bigl(\|q - p_i\|^2\bigr)\phi(q)\,dq = \int_{V_i(P)} 2(p_i - q)\phi(q)\,dq. \]

Split and factor:

\[ \frac{\partial H}{\partial p_i} = 2\left( p_i \int_{V_i(P)} \phi(q)\,dq - \int_{V_i(P)} q\phi(q)\,dq \right) = 2 m_i(P)\bigl(p_i - c_i(P)\bigr), \]

where \( m_i(P) \) and \( c_i(P) \) are the mass and centroid of \( V_i(P) \).

Problem 4 (Discrete-Time Consensus Stability): Consider discrete-time consensus \( z(k+1) = (I_{Nd} - \alpha(L\otimes I_d))z(k) \) with \( \alpha > 0 \). Show that the system is stable (eigenvalues in the closed unit disc, with one eigenvalue at 1) if \( 0 < \alpha < \frac{2}{\lambda_{\max}(L)} \).

Solution:

Since \( L \) is symmetric, it has an eigenbasis \( v_k \) with eigenvalues \( \lambda_k \ge 0 \). The eigenvalues of \( I - \alpha L \) are \( 1 - \alpha\lambda_k \). Stability requires \( |1 - \alpha\lambda_k| \le 1 \) for all \( k \), i.e.

\[ -1 \le 1 - \alpha\lambda_k \le 1 \quad \Rightarrow \quad 0 \le \alpha\lambda_k \le 2. \]

For the largest eigenvalue \( \lambda_{\max}(L) \), this yields \( \alpha \le 2/\lambda_{\max}(L) \), while \( \alpha > 0 \) ensures the nonzero modes decay. The eigenvalue corresponding to \( \lambda_1 = 0 \) is \( 1 \), representing the consensus subspace.

Problem 5 (Desired Distance in Pairwise Potential): For the pairwise potential \( \Phi(r) = \frac{1}{4}(r^2 - d_{\mathrm{des}}^2)^2 \), show that \( r = d_{\mathrm{des}} \) is a local minimum of \( \Phi \) and compute \( \Phi''(d_{\mathrm{des}}) \).

Solution:

We have

\[ \Phi'(r) = (r^2 - d_{\mathrm{des}}^2)r, \]

so critical points occur when \( r = 0 \) or \( r = d_{\mathrm{des}} \). The second derivative is

\[ \Phi''(r) = \frac{d}{dr}\bigl((r^2 - d_{\mathrm{des}}^2)r\bigr) = 3r^2 - d_{\mathrm{des}}^2. \]

Evaluating at \( r = d_{\mathrm{des}} \):

\[ \Phi''(d_{\mathrm{des}}) = 3d_{\mathrm{des}}^2 - d_{\mathrm{des}}^2 = 2 d_{\mathrm{des}}^2 > 0, \]

so \( r = d_{\mathrm{des}} \) is indeed a local minimum, consistent with dispersion around the desired distance.

We developed a unified view of swarm aggregation, dispersion, and coverage from local interaction rules on graphs. Aggregation arises naturally from consensus dynamics driven by graph Laplacians, with Lyapunov functions guaranteeing convergence to the average position. Dispersion is realized via pairwise potentials whose minima encode desired inter-robot distances, leading to shape formation and collision avoidance.

Coverage control relies on Voronoi partitions and quadratic sensing functionals, where gradient descent drives robots toward centroids of their cells, yielding centroidal Voronoi tessellations. Discrete-time stability constraints on step size link continuous-time designs to implementable algorithms. Multilingual implementations (Python, C++, Java, MATLAB/Simulink, Mathematica) illustrate how these theoretical constructs translate into simulation and eventually into controllers for real swarms.