In swarm robotics, we consider a population of \( N \) relatively simple agents, each with a state \( x_i(k) \in \mathbb{R}^d \) at discrete time \( k \in \mathbb{N} \). The collection of all states is

\[ \mathbf{x}(k) = \begin{bmatrix} x_1(k) \\ \vdots \\ x_N(k) \end{bmatrix} \in \mathbb{R}^{Nd}. \]

A local rule specifies how agent \( i \) updates its state based only on information from agents in a local neighborhood

\[ \mathcal{N}_i(k) = \{ j \in \{1,\dots,N\} \mid \|x_j(k)-x_i(k)\| \le R \}, \]

where \( R > 0 \) is an interaction radius. A general local update law is

\[ x_i(k+1) = F_i\!\bigl(x_i(k), \{x_j(k)\}_{j \in \mathcal{N}_i(k)}\bigr), \quad i = 1,\dots,N. \]

An emergent behavior is a global pattern of \( \mathbf{x}(k) \) (e.g. all states converging to a common value, or forming a lattice) that is not explicitly encoded as a centralized command, but arises as the collective effect of all local rules.

Conceptually, we can think of a mapping from micro-rules to macro-patterns: local interaction functions \( F_i \) induce a global dynamical system on \( \mathbf{x}(k) \), whose attractors encode emergent collective behaviors.

flowchart TD
  A["Define local rule for agent i"] --> B["Each agent reads neighbors (within R)"]
  B --> C["Each agent computes control u_i(k)"]
  C --> D["All agents update state in parallel"]
  D --> E["Global swarm state x(k)"]
  E --> F["Observe emergent pattern (e.g. consensus, flocking, formation)"]
        

It is convenient to represent local interactions using a graph \( \mathcal{G}(k) = (\mathcal{V}, \mathcal{E}(k)) \), where \( \mathcal{V} = \{1,\dots,N\} \) is the agent set and \( (j,i) \in \mathcal{E}(k) \) if agent \( j \) sends information to agent \( i \) at time \( k \). The corresponding adjacency matrix \( A(k) = [a_{ij}(k)] \) is defined by

\[ a_{ij}(k) = \begin{cases} 1 & \text{if } (j,i) \in \mathcal{E}(k),\\ 0 & \text{otherwise.} \end{cases} \]

A fundamental class of local rules is averaging dynamics, where each agent moves towards the average of its neighbors:

\[ x_i(k+1) = x_i(k) + \alpha \sum_{j \in \mathcal{N}_i(k)} w_{ij}(k)\bigl(x_j(k)-x_i(k)\bigr), \]

with step size \( \alpha > 0 \) and nonnegative weights \( w_{ij}(k) \ge 0 \). When the graph is fixed and we choose \( W = [w_{ij}] \) such that

\[ w_{ij} \ge 0,\quad \sum_{j=1}^N w_{ij} = 1 \quad\text{for all } i, \]

the stacked dynamics can be written as a linear system

\[ \mathbf{x}(k+1) = (W \otimes I_d)\,\mathbf{x}(k), \]

where \( \otimes \) denotes the Kronecker product and \( I_d \) is the \( d\times d \) identity matrix. The matrix \( W \) is row-stochastic (each row sums to 1) and its graph structure coincides with \( \mathcal{G} \).

A related continuous-time model uses the graph Laplacian \( L = D - A \), where \( D \) is the degree matrix. The dynamics

\[ \dot{\mathbf{x}}(t) = - (L \otimes I_d)\,\mathbf{x}(t) \]

describe consensus flows, which we will interpret as emergent agreement arising from purely local diffusive coupling.

Consensus is a prototypical emergent behavior: all agents converge to a common state even though they only average local information. Consider the discrete-time averaging model with fixed weights:

\[ x_i(k+1) = \sum_{j=1}^N w_{ij} x_j(k), \quad i=1,\dots,N, \]

or in stacked form \( \mathbf{x}(k+1) = (W \otimes I_d)\,\mathbf{x}(k) \). Define the global average

\[ \bar{x}(k) = \frac{1}{N} \sum_{i=1}^N x_i(k). \]

If \( W \) is row-stochastic and \( \mathbf{1} = [1,\dots,1]^\top \), then

\[ \mathbf{1}^\top W = \mathbf{1}^\top, \]

which implies invariance of the average:

\[ \begin{aligned} \bar{x}(k+1) &= \frac{1}{N}\mathbf{1}^\top \mathbf{x}(k+1) = \frac{1}{N}\mathbf{1}^\top (W \otimes I_d)\,\mathbf{x}(k) \\ &= \frac{1}{N} (\mathbf{1}^\top W \otimes I_d)\,\mathbf{x}(k) = \frac{1}{N} (\mathbf{1}^\top \otimes I_d)\,\mathbf{x}(k) = \bar{x}(k). \end{aligned} \]

Thus, the average is a conserved quantity under local averaging.

Suppose further that \( W \) is primitive (it is irreducible and has a strictly positive diagonal). By the Perron–Frobenius theorem, the eigenvalue \( \lambda_1 = 1 \) is simple and all other eigenvalues satisfy \( |\lambda_i| < 1 \). Then

\[ W^k \rightarrow \mathbf{1}\,\pi^\top \quad\text{as } k \rightarrow \infty, \]

where \( \pi \) is the unique stationary distribution of \( W \). Consequently,

\[ \mathbf{x}(k) = (W^k \otimes I_d)\,\mathbf{x}(0) \rightarrow (\mathbf{1}\,\pi^\top \otimes I_d)\,\mathbf{x}(0), \]

which means all agents converge to a common value \( x^\star = \sum_{i=1}^N \pi_i x_i(0) \). The consensus value is an emergent global quantity determined by the initial conditions and the stationary distribution, even though each agent only performs local averaging.

In the special case of a symmetric doubly-stochastic matrix \( W \), we have \( \pi_i = 1/N \), so the emergent consensus equals the initial arithmetic average \( \bar{x}(0) \).

Consensus is linear and leads to all agents collapsing to the same state. Many swarm behaviors require structured configurations: agents should aggregate but maintain spacing, or orbit around a point, etc. A standard approach uses pairwise interaction potentials.

Consider agents in the plane (\( d=2 \)) with positions \( x_i \in \mathbb{R}^2 \). Let a scalar pairwise potential \( \phi : [0,\infty) \rightarrow \mathbb{R} \) depend on distance \( r_{ij} = \|x_i - x_j\| \). A typical attraction–repulsion potential is

\[ \phi(r) = a r^2 - b \log(r+\varepsilon), \quad a,b > 0,\; \varepsilon > 0. \]

The total potential energy of the swarm is

\[ U(\mathbf{x}) = \sum_{1 \le i < j \le N} \phi(\|x_i - x_j\|). \]

A gradient-flow local rule is

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

The interaction is fully local: agent \( i \) uses only the relative positions \( x_j - x_i \) for neighbors \( j \in \mathcal{N}_i \). Yet the global behavior is gradient descent on \( U \):

\[ \frac{\mathrm{d}}{\mathrm{d}t} U(\mathbf{x}(t)) = \sum_{i=1}^N \nabla_{x_i} U(\mathbf{x})^\top \dot{x}_i = - \sum_{i=1}^N \|\nabla_{x_i} U(\mathbf{x})\|^2 \le 0. \]

Thus \( U(\mathbf{x}(t)) \) is non-increasing, and trajectories tend to the set of critical points of \( U \). These critical points correspond to emergent spatial formations (e.g. roughly evenly spaced aggregates), determined by minima of the interaction potential.

flowchart TD
  S["Initial scattered positions x_i(0)"] --> L1["Local distance-based forces from potential phi(r)"]
  L1 --> L2["Agents follow gradient descent of total U(x)"]
  L2 --> L3["Energy U(x(t)) decreases over time"]
  L3 --> L4["Swarm converges to formation (local minimum of U)"]
        

We now implement a simple 2D swarm with local attraction–repulsion using multiple programming languages. The dynamics discretize the gradient-flow rule:

\[ x_i(k+1) = x_i(k) - h \sum_{j \in \mathcal{N}_i(k)} \phi'(\|x_i(k) - x_j(k)\|) \frac{x_i(k) - x_j(k)}{\|x_i(k) - x_j(k)\|}, \]

with step size \( h > 0 \) and potential \( \phi(r) = a r^2 - b \log(r+\varepsilon) \). Neighborhood \( \mathcal{N}_i(k) \) is defined by a sensing radius \( R \).

5.1 Python Implementation (NumPy, basic visualization)

import numpy as np

# Parameters
N = 20          # number of agents
d = 2           # dimension
a = 0.5
b = 1.0
eps = 1e-2
R = 1.0         # sensing radius
h = 0.02        # time step
steps = 1000

# Initialize positions randomly in a square
np.random.seed(0)
x = np.random.uniform(low=-2.0, high=2.0, size=(N, d))

def phi_prime(r):
    # derivative of phi(r) = a r^2 - b log(r + eps)
    return 2.0 * a * r - b / (r + eps)

for k in range(steps):
    x_new = x.copy()
    for i in range(N):
        force = np.zeros(d)
        for j in range(N):
            if i == j:
                continue
            diff = x[i] - x[j]
            dist = np.linalg.norm(diff)
            if dist < 1e-6 or dist > R:
                continue
            fmag = phi_prime(dist)
            force += fmag * diff / dist
        x_new[i] = x[i] - h * force
    x = x_new

# Optionally visualize (requires matplotlib)
if __name__ == "__main__":
    import matplotlib.pyplot as plt
    plt.scatter(x[:,0], x[:,1])
    plt.title("Final swarm configuration")
    plt.axis("equal")
    plt.show()
      

In a robotics stack, this core update could be integrated into a ROS node (e.g., a rclpy or rospy node) that publishes velocity commands to each robot based on locally sensed neighbor positions.

5.2 C++ Implementation (Eigen for Linear Algebra)

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>

using Vec2 = Eigen::Vector2d;

double phi_prime(double r, double a, double b, double eps) {
    return 2.0 * a * r - b / (r + eps);
}

int main() {
    const int N = 20;
    const double a = 0.5;
    const double b = 1.0;
    const double eps = 1e-2;
    const double R = 1.0;
    const double h = 0.02;
    const int steps = 1000;

    std::vector<Vec2> x(N);
    std::vector<Vec2> x_new(N);

    // Simple initialization
    for (int i = 0; i < N; ++i) {
        x = -2.0 + 4.0 * (double)std::rand() / RAND_MAX;
        x = -2.0 + 4.0 * (double)std::rand() / RAND_MAX;
    }

    for (int k = 0; k < steps; ++k) {
        for (int i = 0; i < N; ++i) {
            Vec2 force(0.0, 0.0);
            for (int j = 0; j < N; ++j) {
                if (i == j) continue;
                Vec2 diff = x[i] - x[j];
                double dist = diff.norm();
                if (dist < 1e-6 || dist > R) continue;
                double fmag = phi_prime(dist, a, b, eps);
                force += (fmag / dist) * diff;
            }
            x_new[i] = x[i] - h * force;
        }
        x = x_new;
    }

    // Output final positions
    for (int i = 0; i < N; ++i) {
        std::cout << i << " : " << x[i].transpose() << std::endl;
    }
    return 0;
}
      

This core loop can be embedded in a robotics middleware (e.g., ROS 2 with rclcpp), where positions are updated from sensor messages and control commands are sent via velocity topics.

5.3 Java Implementation (Simple 2D Swarm)

public class LocalRuleSwarm {

    static class Vec2 {
        double x, y;
        Vec2(double x, double y) { this.x = x; this.y = y; }
        Vec2 add(Vec2 other) { return new Vec2(x + other.x, y + other.y); }
        Vec2 sub(Vec2 other) { return new Vec2(x - other.x, y - other.y); }
        Vec2 scale(double s) { return new Vec2(s * x, s * y); }
        double norm() { return Math.sqrt(x * x + y * y); }
    }

    static double phiPrime(double r, double a, double b, double eps) {
        return 2.0 * a * r - b / (r + eps);
    }

    public static void main(String[] args) {
        int N = 20;
        double a = 0.5, b = 1.0, eps = 1e-2;
        double R = 1.0, h = 0.02;
        int steps = 1000;

        Vec2[] x = new Vec2[N];
        Vec2[] xNew = new Vec2[N];

        java.util.Random rng = new java.util.Random(0);
        for (int i = 0; i < N; i++) {
            double px = -2.0 + 4.0 * rng.nextDouble();
            double py = -2.0 + 4.0 * rng.nextDouble();
            x[i] = new Vec2(px, py);
        }

        for (int k = 0; k < steps; k++) {
            for (int i = 0; i < N; i++) {
                Vec2 force = new Vec2(0.0, 0.0);
                for (int j = 0; j < N; j++) {
                    if (i == j) continue;
                    Vec2 diff = x[i].sub(x[j]);
                    double dist = diff.norm();
                    if (dist < 1e-6 || dist > R) continue;
                    double fmag = phiPrime(dist, a, b, eps);
                    Vec2 contrib = diff.scale(fmag / dist);
                    force = force.add(contrib);
                }
                xNew[i] = x[i].sub(force.scale(h));
            }
            System.arraycopy(xNew, 0, x, 0, N);
        }

        for (int i = 0; i < N; i++) {
            System.out.printf("%d : (%.3f, %.3f)%n", i, x[i].x, x[i].y);
        }
    }
}
      

5.4 MATLAB / Simulink Implementation

In MATLAB, we can implement the discrete-time dynamics directly. A Simulink model can use a vector state block with a MATLAB Function block implementing the local rule.

% Parameters
N   = 20;
d   = 2;
a   = 0.5;
b   = 1.0;
eps = 1e-2;
R   = 1.0;
h   = 0.02;
steps = 1000;

% Initialize positions
rng(0);
x = -2 + 4 * rand(N, d);

phi_prime = @(r) 2*a*r - b./(r + eps);

for k = 1:steps
    x_new = x;
    for i = 1:N
        force = zeros(1, d);
        for j = 1:N
            if i == j, continue; end
            diff = x(i,:) - x(j,:);
            dist = norm(diff);
            if dist < 1e-6 || dist > R, continue; end
            fmag = phi_prime(dist);
            force = force + (fmag / dist) * diff;
        end
        x_new(i,:) = x(i,:) - h * force;
    end
    x = x_new;
end

scatter(x(:,1), x(:,2), 'filled');
axis equal;
title('Final swarm configuration');
      

In Simulink, one can create a model with:

• A Discrete-Time Integrator (or Unit Delay plus Sum) holding the stacked positions.
• A MATLAB Function block that takes the current positions x and outputs the vector of increments computed using the same local rule as above.
• Optionally, a visualization subsystem that publishes positions to a robotics visualization tool (e.g., using Robotics System Toolbox for interfacing with external simulation).

5.5 Wolfram Mathematica Implementation

ClearAll["Global`*"];

(* Parameters *)
n = 20;
d = 2;
a = 0.5;
b = 1.0;
eps = 1.*10^-2;
R = 1.0;
h = 0.02;
steps = 1000;

phiPrime[r_] := 2 a r - b/(r + eps);

(* Initialize positions uniformly in a square *)
SeedRandom[0];
x0 = RandomReal[{-2, 2}, {n, d}];

updateSwarm[pts_] := Module[{x = pts, xNew = pts, diff, dist, force},
  Do[
    force = ConstantArray[0., d];
    Do[
      If[i == j, Continue[]];
      diff = x[[i]] - x[[j]];
      dist = Norm[diff];
      If[dist < 10^-6 || dist > R, Continue[]];
      force = force + (phiPrime[dist]/dist) diff;
      ,
      {j, 1, n}
    ];
    xNew[[i]] = x[[i]] - h force;
    ,
    {i, 1, n}
  ];
  xNew
];

x = Nest[updateSwarm, x0, steps];

ListPlot[x,
  PlotStyle -> PointSize[Medium],
  AspectRatio -> 1,
  AxesLabel -> {"x1", "x2"},
  PlotRange -> All,
  PlotLabel -> "Final swarm configuration"
]
      

Problem 1 (Average Invariance in Linear Consensus). Consider the discrete-time consensus dynamics \( x_i(k+1) = \sum_{j=1}^N w_{ij} x_j(k) \) with a fixed row-stochastic matrix \( W \). Show that the arithmetic average \( \bar{x}(k) = \frac{1}{N}\sum_{i=1}^N x_i(k) \) is invariant.

Solution. Stack states as \( \mathbf{x}(k) \) and define \( \mathbf{1} = [1,\dots,1]^\top \). Then

\[ \bar{x}(k) = \frac{1}{N}\mathbf{1}^\top \mathbf{x}(k). \]

Using \( \mathbf{x}(k+1) = W \mathbf{x}(k) \) and row-stochasticity \( \mathbf{1}^\top W = \mathbf{1}^\top \), we get

\[ \bar{x}(k+1) = \frac{1}{N}\mathbf{1}^\top \mathbf{x}(k+1) = \frac{1}{N}\mathbf{1}^\top W \mathbf{x}(k) = \frac{1}{N}\mathbf{1}^\top \mathbf{x}(k) = \bar{x}(k). \]

Hence the average is invariant under the consensus dynamics.

Problem 2 (Lyapunov Function for Consensus). Consider continuous-time consensus \( \dot{\mathbf{x}} = - (L \otimes I_d)\,\mathbf{x} \) on a connected undirected graph with Laplacian \( L \). Define

\[ V(\mathbf{x}) = \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N a_{ij} \|x_i - x_j\|^2. \]

Show that \( V(\mathbf{x}(t)) \) is non-increasing along trajectories.

Solution. First note that for an undirected graph \( a_{ij} = a_{ji} \). Differentiating,

\[ \frac{\mathrm{d}}{\mathrm{d}t} V = \sum_{i=1}^N \sum_{j=1}^N a_{ij} (x_i - x_j)^\top (\dot{x}_i - \dot{x}_j). \]

Using \( \dot{x}_i = -\sum_{\ell=1}^N a_{i\ell}(x_i - x_\ell) \), we get after algebraic simplification the identity

\[ \frac{\mathrm{d}}{\mathrm{d}t} V = - \sum_{i=1}^N \left\|\sum_{j=1}^N a_{ij} (x_i - x_j) \right\|^2 \le 0. \]

Thus \( V \) is a Lyapunov function, non-increasing along trajectories. Since \( V \ge 0 \), it converges and the swarm approaches the largest invariant set where the right-hand side vanishes, which corresponds to \( x_i = x_j \) for all edges, i.e. consensus.

Problem 3 (Equilibrium Distance in Attraction–Repulsion). For the potential \( \phi(r) = a r^2 - b \log(r+\varepsilon) \) with \( a,b > 0 \), compute the positive equilibrium distance \( r^\star \) between two agents, i.e. the distance where \( \phi'(r^\star) = 0 \).

Solution. Compute

\[ \phi'(r) = 2 a r - \frac{b}{r+\varepsilon}. \]

Set \( \phi'(r^\star) = 0 \):

\[ 2 a r^\star - \frac{b}{r^\star+\varepsilon} = 0 \quad \Rightarrow \quad 2 a r^\star (r^\star+\varepsilon) = b. \]

This yields a quadratic equation

\[ 2 a (r^\star)^2 + 2 a \varepsilon r^\star - b = 0. \]

Solving, we obtain

\[ r^\star = \frac{-2 a \varepsilon + \sqrt{4 a^2 \varepsilon^2 + 8 a b}}{4 a} = -\frac{\varepsilon}{2} + \frac{1}{2} \sqrt{\varepsilon^2 + \frac{2 b}{a}}. \]

The other root is negative and thus discarded. The swarm tends to configurations where typical inter-agent distances are near \( r^\star \).

Problem 4 (Connectivity and Consensus). Let \( W \) be a row-stochastic matrix with graph \( \mathcal{G} \). Assume that \( \mathcal{G} \) is strongly connected and \( w_{ii} > 0 \) for all \( i \). Explain (without a full proof) why the consensus value is unique.

Solution. Strong connectivity and positive diagonal ensure that \( W \) is primitive. By Perron–Frobenius, 1 is a simple dominant eigenvalue, and there exists a unique positive left eigenvector \( \pi^\top W = \pi^\top \) with \( \sum_i \pi_i = 1 \). For any initial condition, \( W^k \rightarrow \mathbf{1}\,\pi^\top \), so

\[ \mathbf{x}(k) \rightarrow (\mathbf{1}\,\pi^\top) \mathbf{x}(0). \]

The asymptotic consensus value is thus uniquely determined as \( x^\star = \sum_{i=1}^N \pi_i x_i(0) \), independent of time and the specific trajectory, depending only on \( \pi \) and the initial condition.

Problem 5 (Energy Descent for Discrete-Time Potential Dynamics). Consider the discrete-time update \( x_i(k+1) = x_i(k) - h \nabla_{x_i} U(\mathbf{x}(k)) \) with small \( h > 0 \), where \( U(\mathbf{x}) \) is a smooth potential. Show that for sufficiently small \( h \), the potential \( U(\mathbf{x}(k)) \) decreases along iterations, up to first order in \( h \).

Solution. Using first-order Taylor expansion,

\[ U(\mathbf{x}(k+1)) \approx U(\mathbf{x}(k)) + \nabla U(\mathbf{x}(k))^\top (\mathbf{x}(k+1) - \mathbf{x}(k)). \]

But \( \mathbf{x}(k+1) - \mathbf{x}(k) = - h \nabla U(\mathbf{x}(k)) \). Thus

\[ U(\mathbf{x}(k+1)) \approx U(\mathbf{x}(k)) - h \|\nabla U(\mathbf{x}(k))\|^2. \]

Hence, for small \( h \), the potential decreases by approximately \( h \|\nabla U(\mathbf{x}(k))\|^2 \) at each step, unless the gradient vanishes. More precise analysis uses higher-order terms to bound the step size that guarantees monotone decrease.

In this lesson we modeled swarm robotics as a networked dynamical system in which each agent obeys simple local rules based on neighbor information. We showed how linear averaging dynamics yield emergent consensus, with invariance of the global average and convergence characterized by Perron–Frobenius theory. We then introduced nonlinear interaction potentials, where gradient-flow local rules cause the swarm to minimize a global energy function and organize into spatially structured formations. Finally, we implemented these ideas in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, providing a computational bridge from theory to simulation.