Distributed Predictive Flocking for a Swarm of Miniature Vibration-Driven Robots in a Confined Arena

_gset:Nn __short_title: Distributed Predictive Flocking for a Swarm of Miniature Vibration-Driven Robots in a Confined Arena

_gset:Nn __short_authors: Shrestha et al.

[orcid=0009-0005-9585-3800]

Introduction

Autonomous multi-robot systems are a promising approach for tasks in cluttered, inaccessible, or hazardous spaces, from confined-structure inspection to search in collapsed buildings and medical micro-robotics [14]. In these settings, many miniature robots must coordinate in proximity within a tight, bounded arena. As in natural swarms, such as flocking of starlings, schooling of fish, the crowded colonies of ants and honey bees, and even the collective migration of cells, robust group behavior emerges from purely local interactions among simple agents [9, 7, 6], with no central control [14]. Miniaturization, however, tightly bounds onboard compute, memory, and power [15], demanding control that is theoretically principled yet computationally light. Moreover, when the arena is fixed, a growing number of robots increases the density of the swarm and hinders smooth motion. When it is also cluttered with obstacles, coordination is even harder. In this scenario, collision avoidance becomes the central challenge.

The standard way to coordinate such a swarm from purely local interactions is Reynolds’ boid model [1, 2], which reduces collective motion to three local rules (cohesion, alignment, and separation) that produce coherent group motion without global information [8]. Early control laws formalized these as artificial potential fields over a communication graph [4, 5], and reactive velocity-matching rules scale to tens of drones [13]. Such reactive methods, however, offer no predictive horizon, no formal safety guarantee under nonlinear dynamics, and no systematic handling of kinodynamic limits. These shortcomings are especially severe in confined, high-density settings, where many underactuated robots have little room to recover from unsafe motions.

Model predictive control (MPC) [18] addresses these limitations. For the nonlinear dynamics of mobile robots, this takes the form of nonlinear MPC (NMPC). It solves a finite-horizon optimal control problem (OCP) online, which provides closed-loop feedback, constraint satisfaction, and explicit multi-objective trade-offs. Distributed MPC (DMPC) extends this to a swarm by splitting the global problem into per-agent OCPs coupled through shared constraints [19, 20, 21, 22]. Solvers such as the Sequential Quadratic Programming real-time iteration (SQP-RTI) [55] solve each per-agent OCP in real time. Similarly, optimization-based swarm planners achieve strong real-world navigation [23, 25, 24, 27]. These methods are mainly developed for point-to-point navigation [20, 26], with collision avoidance typically encoded as geometric constraints such as separating hyperplanes or safe corridors.

The closest works embed Reynolds behaviors in NMPC. They have been demonstrated on up to 16 quadrotors [28, 29], leader-follower teams [30, 31, 32], and fixed-wing unmanned aerial vehicles (UAVs) [33]. Most share one regime. They assume sparse flight in large 3D airspace, with the swarm migrating toward a shared goal. Even when several obstacles are present, they remain sparse relative to the available space, keeping the environment effectively low-density. Few address the opposite regime that we target, where many underactuated robots coordinate at high density in a fixed, confined, and cluttered arena under a task objective such as following a light gradient. Here, the inter-robot safety formulation, not free-space planning, governs how many robots can coexist. Each of these works, however, fixes a single safety formulation in advance, usually a distance constraint or a penalty. It remains open which formulation fits the most robots into a fixed arena with minimal collisions.

Safety formulations differ in how collision avoidance is encoded as a constraint in the OCP. Such a constraint is classically formulated using a barrier function $h(x)$ whose safe set $\mathcal{S} = \{x \mid h(x)\ge 0\}$ is kept forward invariant to certify safety [34], typically imposed as a hard [36] or soft [37] distance constraint. A Control Barrier Function (CBF) [34] relaxes this strict boundary condition by allowing $h(x)$ to decrease at a rate bounded by a class-$\mathcal{K}$ function, and its discrete-time form (Dt-CBF) [35] imposes this at each MPC stage. [38] enforces it as a hard constraint, which is often infeasible. Feasibility is restored by softening the decay with slack variables [39, 40, 41], adding terminal barrier terms [43, 44], or anchoring a soft CBF with one hard step [42]. Beyond feasibility, when the constraint has relative degree greater than one, as for a distance constraint under unicycle dynamics, a higher-order CBF (HOCBF) [45, 46] chains the barrier conditions. [42] compare a distance-constraint MPC, a hard discrete-time CBF, a slacked CBF, and their combination on a common footing, but for a single robot. In the multi-robot setting, safety is enforced by reactive schemes [48, 49], by post-hoc CBF-QP (quadratic program) filters [50], or by a CBF inside a decentralized MPC with neighbors as predicted trajectories [51, 52, 53, 54]. A formulation’s collision performance, however, likely changes with the swarm density. A comparison across densities is required to find the best one.

In this work, we investigate the maximum number of miniature robots a distributed predictive controller can flock in a fixed arena, and design a controller to push that limit. Each robot solves its own NMPC combining Reynolds cohesion, alignment, and separation (gated by neighbor presence) with inter-robot/obstacle safety and the task objective. It selects the $K$ nearest neighbors within the communication radius and exchanges predicted trajectories synchronously as OCP parameters, giving a fully distributed, leader-free architecture with no central planner or global map. We evaluate it on a physics-informed digital twin [17] of the GRASPion [16], an open-source vibration-driven bristlebot studied as an active-matter model system. We characterize how densely the swarm packs a $0.64\,\mathrm{m^2}$ arena while sustaining two behaviors. These are phototaxis (ascending toward a moving Gaussian light gradient) and a phototactic orbital mode (orbiting the source along the intensity level-set). We compare six inter-robot safety formulations on a common footing (§ 3.3). A slacked barrier paired with a discrete-time CBF at a forward handle [47] sustains the highest density with the fewest collisions. The main contributions of this paper are as follows:

1. A distributed predictive control system for dense swarms of miniature vibration-driven robots. Each robot solves its own NMPC fusing flocking, safety, and the task objective, coupled to neighbors only through exchanged predicted trajectories.

2. Demonstration that a single controller generalizes across qualitatively different collective behaviors. The same architecture produces gradient-ascending phototaxis, stable phototactic orbital formation, and generic flocking with no modification to the safety or coordination layers, showing that task objectives and coordination structure are separable and that new tasks can be added without re-deriving safety guarantees.

3. A careful comparison of six inter-robot safety formulations within one identical controller, isolating the safety constraint as the only varying factor. The Look-ahead Slack Dt-CBF gives the best overall trade-off for the phototaxis deployment target.

4. Characterization of the operating swarm density that the controller sustains in a fixed arena. The swarm stays cohesive with few collisions in small swarms; collisions rise steeply towards the geometric packing limit, as saturation sets in, while coordination and feasibility degrade gently rather than collapse.

The remainder of this paper is structured as follows. § 2 presents the robot model, light field, control architecture overview, behavior generation, flocking rules, and the DMPC OCP. § 3 reports the simulation results, including the safety-formulation study (§ 3.3), the analysis of collective behavior (§ 3.4), and the quantification of collision avoidance under different tasks (§ 3.5). § 4 summarizes the work and provides an outlook for future studies.

Preliminaries and Problem Statement

In this work, we consider a swarm of $N_b$ robots operating in a bounded arena $\mathcal{A} = [-A_h,\,A_h]^2 \subset \mathbb{R}^2$, with each robot labeled by $i \in \{1,\ldots,N_b\}$. The state of robot $i$ is $\bm{x}_i = [\bm{p}_i,\,\theta_i]^\top \in \mathbb{R}^3$, where $\bm{p}_i = [p_{x,i},\,p_{y,i}]^\top$ is the position and $\theta_i$ the heading. Inter-robot communication is modeled as a directed sensing graph $\mathcal{G} = (\mathcal{V},\,\mathcal{E})$, where $\mathcal{V} = \{1,\ldots,N_b\}$ denotes the set of robots and $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$ is the communication set of ordered pairs $(i,j)$ for which robot $i$ can observe robot $j$. The neighbor set of robot $i$ is $\mathcal{N}_i = \{j \in \mathcal{V} \mid (i,j) \in \mathcal{E}\}$. We adopt a topological neighborhood with fixed cardinality $|\mathcal{N}_i| \le K$, selecting the $K$ nearest robots within communication radius $R_{\rm comm}$ of robot $i$. This choice keeps the neighbor count constant and avoids the cost discontinuities that arise under Voronoi or attraction-based topologies [10, 11, 12]. Each robot additionally accounts for up to $M$ nearby static circular obstacles within a sensing range $d_{\rm obs}$, indexed by $m \in \{1,\ldots, M\}$, with obstacle $m$ at center $\bm{p}_m \in \mathbb{R}^2$ and radius $r_m$. The environment also includes a single light source whose intensity field serves as the swarm stimulus.

The OCP is compiled once with a fixed parameter dimension. Neighbor and obstacle information therefore enters two fixed-size blocks of $K$ and $M$ slots, the neighbor block and the obstacle block. At each control step, robot $i$ receives a predicted trajectory $\bm{\tau}_j = \bigl\{[\hat{\bm{p}}_j^{(k)},\,\hat{\bm{v}}_j^{(k)}]^\top\bigr\}_{k=0}^{N}$ from each neighbor $j \in \mathcal{N}_i$, where $\hat{\bm{p}}_j^{(k)} \in \mathbb{R}^2$ and $\hat{\bm{v}}_j^{(k)} \in \mathbb{R}^2$ are the predicted position and velocity of neighbor $j$ at horizon stage $k$. Throughout, a hat denotes a quantity the per-robot OCP treats as fixed data (the ego robot’s current measured state and the neighbors’ received predicted trajectories). The ego robot’s own predicted states are decision variables. These trajectories fill the neighbor block, and the number of occupied slots changes from step to step. Each neighbor slot therefore carries a binary activation mask $a_j^{(k)} \in \{0,1\}$, where $a_j^{(k)} = 1$ if slot $j$ is occupied at stage $k$ and $a_j^{(k)} = 0$ otherwise. The obstacle block is filled the same way. Robot $i$ selects the $M$ nearest obstacles within obstacle sensing range $d_{\rm obs}$ and writes each center $\bm{p}_m$ into a slot together with its own activation flag.

Robot Model

Each robot is modeled as a rigid ellipse with semi-major axis $a$ and semi-minor axis $b$, matching the physical footprint of the GRASPion open-source bristlebot [16], as illustrated in Figure 1. The simulated elliptical agent is modeled after this hardware platform, following the physics-informed digital twin of [17], so the simulation study is grounded in a physical robot footprint and can serve as a precursor to future physical experiments.

The dynamics of each robot are described by a three-state kinematic unicycle model with state $\bm{x} \in \mathbb{R}^3$, input $\bm{u} \in \mathbb{R}^2$, and continuous-time dynamics $f(\bm{x}, \bm{u})$ given by

where $(p_x, p_y)$ is the position, $\theta$ the heading angle, $v$ the linear velocity, $v_x = v\cos\theta$ and $v_y = v\sin\theta$ are the velocity components, and $\omega$ the angular velocity. The continuous-time dynamics (3) are integrated using an explicit Runge-Kutta 4 (ERK-4) integrator with three steps per control interval, yielding the discrete-time map $\bm{x}^{(k+1)} = F(\bm{x}^{(k)}, \bm{u}^{(k)})$.

Light Field Model

The arena contains a time-varying light field that each robot senses locally to compute its navigation reference. The field is modeled as a Gaussian intensity profile whose source either orbits the arena center or remains stationary.

Light Source

The light source produces a Gaussian intensity field

where $\bm{p} \in \mathbb{R}^2$ is an arbitrary position, $I_0 > 0$ is the peak intensity, $\bm{s}(t) \in \mathbb{R}^2$ is the position of the source at time $t$, and $\sigma_{\rm L}$ represents the spatial spread of intensity.

Source Dynamics

The light source supports two motion modes. In circular rotation mode, the source orbits the center of the arena at a fixed radius $r$ with angular speed $\omega_s = v_s /r$, with the position of the source given by

where $\varphi_0$ is the initial angle of the source in the orbit and $v_s$ is the speed of the source along the orbit. In stationary mode, the source remains fixed at the arena center.

Intensity Gradient

The intensity gradient with respect to the robot position is computed analytically as

This analytic gradient is the local cue that the collective-behavior modes of § 2.4 convert into each robot’s look-ahead reference.

Control Architecture Overview

Every robot runs an identical decentralized controller, outlined in Figure 2. At each control step, the robot senses the local light field, and a behavior module maps this reading to a look-ahead reference that encodes where the robot should head next. A per-robot optimal control problem (OCP) tracks that reference. Flocking rules keep the group cohesive and heading-aligned, while barrier and control-barrier-function constraints enforce inter-robot and obstacle safety. Robots publish their predicted trajectories over the sensing graph and re-solve synchronously from their neighbors’ previous plans, coordinating the swarm without a central node. The rest of this section develops the remaining components in turn. These are the behavior module (§ 2.4), the flocking rules (§ 2.5), the six safety formulations (§ 2.5.3–§ 2.5.4), the reference-tracking and control-effort costs (§ 2.6), the assembled OCP (§ 2.7), and the decentralized execution loop (§ 2.8).

Collective Behavior Modes

Collective motion emerges when individual agents produce coordinated group-level movement through purely local interactions in the absence of centralized planning [1, 3]. In our controller, a behavior module translates local sensory observations into a look-ahead reference $(\bm{p}^{\rm ref}, \theta^{\rm ref}, v^{\rm ref}, \omega^{\rm ref})$ at each control step, which the OCP then tracks while maintaining flocking coherence and safety. We describe two such modes, phototaxis and phototactic orbital behavior.

At each control step, the light-field gradient $\bm{g}=\nabla I(\hat{\bm{p}}_i)$ is evaluated analytically from (4)–(6). The resulting reference heading depends on the selected collective behavior mode. Once the heading is obtained, the look-ahead reference position is computed as

where $d_{\rm la}$ is the look-ahead distance, and the clipping function constrains the reference position to remain inside the arena. The reference velocity is set to a constant target value $v^{\rm ref}$, and the reference angular velocity $\omega^{\rm ref}$ is fixed at zero.

Phototaxis Behavior Mode

Phototaxis drives each robot to ascend the intensity gradient of the light field without knowledge of the global source position. The reference heading is aligned with the steepest ascent direction:

where $\epsilon_{\rm grad}$ is a small gradient threshold. When the field is nearly flat ($\|\bm{g}\|_2 < \epsilon_{\rm grad}$), the robot holds its current heading $\hat{\theta}_i$ rather than turning randomly.

Phototactic Orbital Behavior Mode

In the phototactic orbital mode, each robot follows the tangent of the light intensity level set, effectively orbiting the source without knowledge of its global position. The reference heading is obtained by rotating the gradient direction from (8) by $\pi/2$, with the orbital direction governed by the circulation parameter $\kappa \in \{-1, +1\}$ (where $\kappa = +1$ generates counter-clockwise motion, and $\kappa = -1$ generates clockwise motion):

Extensibility to Arbitrary Behaviors

Decoupling the reference generator from the OCP makes the architecture task-agnostic. Any behavior (directed migration, dispersal, transport, shape formation) can be integrated by mapping local observations to a look-ahead reference $(\bm{p}^{\rm ref}, \theta^{\rm ref}, v^{\rm ref}, \omega^{\rm ref})$, while the DMPC preserves safety and flocking coherence.

Reynolds Flocking Rules

The behavior modes set each robot’s individual objective, while the collective motion arises from the local steering rules proposed by Reynolds [1, 2]. Cohesion binds robots into a group, alignment coordinates their heading, separation keeps them from colliding, and obstacle avoidance steers them clear of static hazards. The controller encodes cohesion and alignment as quadratic cost terms in the per-robot OCP, while separation and obstacle avoidance are enforced through barrier and discrete-time CBF constraints [34, 35]. We focus on the separation and obstacle avoidance by comparing six different formulations, ranging from hard barriers to slacked variants (§ 2.5.3 and § 2.5.4).

Cohesion

The active-neighbor centroid at stage $k$ is computed as a masked weighted average of predicted neighbor positions:

where the soft neighbor count

prevents division by zero when there are no neighbors in the vicinity. The cohesion penalty is

with weight $q_{\rm coh} \ge 0$, gated by the activation gate

When at least one neighbor slot is active, $g^{(k)}=1$ and the cost is enabled. In isolation, when $g^{(k)}=0$ and it vanishes. Separating $s^{(k)}$ (denominator safety) from $g^{(k)}$ (which enables the cohesion and alignment costs only when a neighbor is present) keeps the per-robot cost well posed for the SQP-RTI solve. Finally, cohesion uses the stage-matched predicted positions $\hat{\bm{p}}_j^{(k)}$ rather than current measured ones, pulling robot $i$ toward where its neighbors are predicted to be.

Alignment

The mean heading direction ($\theta_{\rm mean}^{(k)}$) of active neighbors is estimated from their predicted velocity vectors:

The neighbor velocities are summed rather than normalized, so faster neighbors weigh more heavily. The $\varepsilon$ offsets regularize $\operatorname{atan2}$ near a zero velocity sum (fixing $\theta_{\rm mean}^{(k)}$ at $\pi/4$ only at exact cancellation, an edge case). The resulting $\theta_{\rm mean}^{(k)}$ then defines the alignment heading error, wrapped to $(-\pi,\pi]$:

and the alignment penalty is

with weight $q_{\rm ali} \ge 0$. As in the cohesion term, the shared gate $g^{(k)}$ disables both cohesion and alignment when robot $i$ is isolated.

Inter-Robot Separation

The inter-robot separation constraint builds on a pairwise safety barrier and its Dt-CBF extension. For neighbor slot $j$ at stage $k$, the barrier

is non-negative exactly when robot $i$ and neighbor $j$ are at least $D_{\rm safe}$ apart. Starting from this barrier, the same safety constraint can be encoded in several ways, each imposed hard or slacked. We consider six formulations that differ only in the safety constraint, and observe the corresponding change in swarm behavior. § 3.3 then evaluates all six under identical conditions and retains the best-performing one. The last three share a common structure consisting of a slacked barrier anchor paired with a rate condition. They differ only in the rate condition, which is a one-step Dt-CBF on the body, the same condition on a forward handle, or a two-step higher-order chain.

(a) Hard BF

The plain barrier imposed as a hard constraint,

is the most direct encoding of the $D_{\rm safe}$ rule. It forbids any predicted configuration closer than $D_{\rm safe}$ at every stage. It is purely reactive (it acts only once a neighbor is within range) and, as a hard constraint, allows no slack, so a robot surrounded by several neighbors at once can leave its OCP with no feasible input.

(b) Hard Dt-CBF

The discrete-time CBF condition as a hard constraint,

compares the barrier $h_j^{(k)}$ with its one-step-ahead value $h_j^{(k+1)}$. The latter is formed by a forward-Euler step of the ego position and a constant-velocity step of each neighbor,

The ego step is linear in $v^{(k)}$, which keeps the quadratic-program subproblems well conditioned. It is a first-order approximation used only inside the barrier constraint, separate from the ERK-4 discretization of the dynamics (3). The condition does not forbid $h_j<0$ but limits how fast the barrier may shrink each step, to at most a fraction $\alpha$ (a discrete class-$\mathcal{K}$ rate). The robot therefore starts slowing while still safely separated rather than braking at the boundary, giving smoother, early avoidance. This is exactly the MPC-CBF formulation of [38]. Being a hard constraint, it carries the same infeasibility risk as Hard BF with dense congregations of robots.

(c) Slack BF

The plain barrier mentioned in Hard BF, softened by its slack,

keeps the OCP solvable in the crowded configurations that make Hard BF infeasible. When no input can satisfy $h_j\ge0$, the penalized slack $\xi_{j,\mathrm{bf}}$ opens just enough to admit a feasible, minimally-violating solution instead of failing.

(d) Slack Dt-CBF

The slacked barrier anchor (23) and the discrete-time CBF condition softened by its own slack,

imposed together on the same neighbor, each with its own slack. The rate condition combines the before-contact activation of Hard Dt-CBF with the improved solvability of Slack BF. The anchor guards the footprint at every stage while the CBF constraint shapes the approach rate. This is the reference point for the next two formulations: they retain this anchored structure and change only where or how the rate condition is evaluated, so the ensuing description only mentions those two choices.

(e) Look-ahead (LA) Slack Dt-CBF

This formulation is a modification of Slack Dt-CBF with its rate condition moved from the body onto a forward handle a fixed distance $L$ ($0<L<D_{\rm safe}$) ahead of the robot, adapting the offset-point construction of [47] to discrete time,

with the slacked condition $h_{j,L}^{(k+1)} - (1-\alpha)\,h_{j,L}^{(k)} + \xi_{j,\mathrm{cbf}}^{(k)} \ge 0$ imposed on the handle. A body distance barrier has relative degree two under unicycle dynamics, so the turn rate does not appear in its one-step form and a robot can only slow down to keep clear. The handle $\bm{p}_L^{(k)}$ moves with the heading, so $\omega^{(k)}$ enters the one-step barrier directly. The robot can, therefore, also steer around a neighbor, while the body anchor (23) still guards the true footprint.

(f) Slack Dt-HOCBF

The slacked barrier anchor (23) together with a slacked higher-order (relative-degree-two) chain [46, 45] is built natively in discrete time from successive barrier values rather than from derivatives. With $h_j^{(k+1)}$ and $h_j^{(k+2)}$, the barrier evaluates the one-step and two-step prediction under the held input $(v^{(k)},\omega^{(k)})$,

with decay rates $\gamma_1,\gamma_2\in(0,1]$. Expanding the chain gives the single constraint actually imposed,

where the two-step prediction advances the heading as $\theta^{(k+1)} = \theta^{(k)} + \Delta t\,\omega^{(k)}$ and each neighbor by its constant velocity. This addresses the same relative-degree-two obstruction as (e) but analytically rather than geometrically: $\omega^{(k)}$ is absent from the one-step barrier, yet it enters $h_j^{(k+2)}$ through the rotated heading, with sensitivity $\partial\psi_2/\partial\omega \propto \Delta t^{2} v^{(k)}$.

Across the six formulations, a safety constraint is either imposed hard or softened by its own non-negative slack ($\xi_{j,\mathrm{bf}}^{(k)}$ for a barrier constraint, $\xi_{j,\mathrm{cbf}}^{(k)}$ for a Dt-CBF constraint), and the active slacks are penalized in the cost,

where $\bm{\xi}_j^{(k)}$ collects the chosen formulation’s inter-robot slacks for neighbor $j$, and $Z^{\rm rob},z^{\rm rob}$ are the quadratic and linear weights. The set $\bm{\xi}_j^{(k)}$ holds both slacks for the hybrid, look-ahead, and higher-order variants (d)–(f), a single slack for Slack BF, and none for the hard variants. Each slack is driven to zero whenever its constraint is satisfiable and absorbs only the unavoidable residual otherwise. Inactive neighbor slots are padded ($h\equiv100$) so their constraints and slacks vanish.

Obstacle Avoidance

The selected formulation (a)–(f) from § 2.5.3 is applied unchanged to static obstacles, so no new safety constraints are introduced. Each obstacle $m$ is treated as a neighbor with a fixed prediction. Three substitutions turn each inter-robot safety constraint of the chosen formulation into its obstacle counterpart. The neighbor position $\hat{\bm{p}}_j^{(k)}$ becomes the obstacle center $\bm{p}_m$ (held at zero velocity). The clearance $D_{\rm safe}$ becomes $D_{\rm safe, obs}$ (a fixed radius exceeding the obstacle radius $r_m$). The class-$\mathcal{K}$ gains take their obstacle values (defaulting to the inter-robot ones). The obstacle barrier, for instance, is

The obstacle slacks are penalized analogously,

with $\bm{\xi}_m^{(k)}$ the chosen formulation’s obstacle slacks for obstacle $m$ and $Z^{\rm obs},z^{\rm obs}$ the quadratic and linear weights. At the terminal stage $k=N$, the one-step prediction is undefined, and all formulations revert to a plain terminal barrier.

Look-Ahead Reference Point and Control Effort

The look-ahead reference point $(\bm{p}^{\rm ref}, \theta^{\rm ref}, v^{\rm ref}, \omega^{\rm ref})$ produced by the behavior module of § 2.4, whatever the mode, enters the per-robot OCP through a quadratic tracking cost. The stage-$k$ deviations of the state from this reference point are then defined as

The wrapping in (34) maps the heading deviation to $(-\pi,\pi]$, avoiding discontinuities at $\pm\pi$. The deviations form a quadratic reference-tracking cost:

with weights $q_{\rm pos}, q_{\theta}, q_v, q_{\omega} > 0$. A control-effort term regularizes the inputs for smoother, more energy-efficient trajectories:

where $r_u > 0$.

OCP Formulation

Each robot solves a finite-horizon OCP over $N$ stages. The cost assembles the look-ahead, control-effort, cohesion, and alignment terms plus CBF slack penalties into a stage and a terminal cost. This cost is minimized subject to the dynamics, the state and input box constraints, and safety constraints that enforce separation and obstacle avoidance.

Stage Cost

The stage cost at horizon stage $k \in \{0,\ldots,N-1\}$ is

where the look-ahead, control-effort, cohesion, alignment, inter-robot slack, and obstacle slack terms are given by (36), (37), (12), (17), (30), and (32), respectively. Separation and obstacle avoidance no longer appear as cost terms. They are enforced through the barrier and Dt-CBF constraints of § 2.5.3–§ 2.5.4.

Terminal Cost

The terminal cost at stage $k = N$ is

where the position, heading, inter-robot slack, and obstacle slack terms are given by (33), (34), (30), and (32), respectively.

Full OCP

Every robot $i$ solves the following optimization problem (with neighbor data treated as stage-wise parameters), which incorporates all the previously defined cost terms and constraints in a unified formulation:

Here, $\bm{C}_{\rm rob\text{-}safety}$ and $\bm{C}_{\rm obs\text{-}safety}$ are the inter-robot and obstacle safety constraints of the selected formulation (a)–(f) (§ 2.5.3 and § 2.5.4), evaluated on the neighbors’ predicted trajectories $\{\bm{\tau}_j\}$ and the obstacle positions $\{\bm{p}_m\}$. The bounds $(\bm{x}_{\min},\bm{x}_{\max})$ and $(\bm{u}_{\min},\bm{u}_{\max})$ are the state and input box limits. The slack vectors $\bm{\xi}_{\rm rob}^{(k)}$ and $\bm{\xi}_{\rm obs}^{(k)}$ hold up to $2K$ and $2M$ entries for $k<N$. The predictive Dt-CBF rows require the one-step-ahead prediction, while the higher-order rows require the two-step prediction, so they drop at the terminal stage $k=N$. For the hard formulations (a) and (b), the slacks are fixed to $\bm{0}$.

Simulation Workflow

Each control step executes a synchronous parallel update [59] in which all robots compute their optimal control sequences simultaneously using only neighbor information from the previous timestep (Algorithm alg:dmpc). This structure prevents any robot from exploiting another robot’s within-step update, matching the limited communication of a real distributed system. At the beginning of each control step, a KD-Tree [60] is constructed over the current robot positions $\{\hat{\bm{p}}_1, \ldots, \hat{\bm{p}}_{N_b}\}$. This spatial partitioning structure enables all $K$-nearest-neighbor queries within communication range $R_{\rm comm}$ to be performed in $O(N_b \log N_b)$ total time compared to the naive pairwise distance search $O(N_b^2)$ [61]. Each robot then samples the local light intensity and gradient to compute a look-ahead reference $(\bm{p}^{\rm ref}, \theta^{\rm ref}, v^{\rm ref}, \omega^{\rm ref})$ which encodes the desired position, heading, linear speed, and angular velocity according to the active swarm behavior mode. It also collects the nearby obstacle positions within range $d_{\rm obs}$ and the neighbors’ predicted trajectories $\{\bm{\tau}_j\}$ buffered from the previous control step. Using these inputs, each robot solves a local OCP over horizon $N$ via a single SQP-RTI step, while updating the stage-wise parameter vector $\bm{\pi}^{(k)}$ to incorporate the look-ahead reference, neighbor trajectories, and obstacle positions as time-varying coupling parameters in (?). After optimization, the robot applies the first control action $\bm{u}^{(0)}$ and publishes its updated predicted trajectory $\bm{\tau}_i$ for use by neighboring robots at time $t + \Delta t$.

\small
  \caption{DMPC Control Step at time $t$}
  \label{alg:dmpc}
  \begin{algorithmic}[1]
    \State Advance light field: $I(\cdot,t) \leftarrow I(\cdot,t{+}\Delta t)$
    \State \textbf{Spatial Indexing:} Construct KD-Tree from current positions
    $\{\hat{\bm{p}}_1, \ldots, \hat{\bm{p}}_{N_b}\}$

    \For{robot $i = 1,\ldots,N_b$ \textbf{in parallel}}
    \Comment{Communication Phase}
    \State Detect $K$ nearest neighbors $j \ne i$ via KD-Tree query with range
    $R_{\rm comm}$
    \State Buffer $\bm{\tau}_j = \{[\hat{\bm{p}}_j^{(k)},
    \hat{\bm{v}}_j^{(k)}]\}_{k=0}^{N}$ from neighbors'
    \emph{previous} trajectories
    \EndFor

    \For{robot $i = 1,\ldots,N_b$ \textbf{in parallel}}
    \Comment{Computation Phase}
    \State Sense light intensity $I(\hat{\bm{x}}_i)$ and gradient
    $\nabla I(\hat{\bm{x}}_i)$
    \State Compute look-ahead reference
    $(\bm{p}^{\rm ref},\, \theta^{\rm ref},\, v^{\rm ref},\, \omega^{\rm ref})$
    based on the active swarm behavior
    \State Detect the $M$ nearest obstacles within range $d_{\rm obs}$ and
    compute their center positions $\bm{p}_m$
    \State Fetch trajectories $\bm{\tau}_j$ from neighbor buffer
    \State Warm-start: $\bm{x}^{(k)} \leftarrow \hat{\bm{x}}_i$,\;
    $\bm{u}^{(k)} \leftarrow \bm{u}_{\rm nom}$ $\forall k$
    \For{$k = 0,\ldots,N$}
    \State Build and set parameter block $\bm{\pi}^{(k)}$: look-ahead reference +
    neighbor trajectory $\bm{\tau}_j$
    + obstacle positions $\bm{p}_m$
    \EndFor
    \State Execute one SQP-RTI step to obtain $\bm{u}^{(k)}$
    \State Extract first optimal control input $\bm{u}_i^* \leftarrow \bm{u}^{(0)}$ and updated
    $\bm{\tau}_i$ from solver
    \EndFor

    \State \textbf{Synchronize:} All robots publish updated trajectories $\bm{\tau}_i$
    simultaneously \Comment{Available at $t{+}\Delta t$}

    \For{robot $i = 1,\ldots,N_b$}
    \Comment{Actuation Phase}
    \State Apply $\bm{u}_i^*$
    \EndFor
  \end{algorithmic}

Simulation Results

We evaluate the proposed DMPC framework through a series of simulation studies. The studies characterize the operating density up to the geometric packing limit, compare the six inter-robot safety formulations under a common controller setup, demonstrate emergent collective behavior across the task modes, and assess real-time computational scalability with swarm size.

Simulation Setup

We implemented our DMPC framework symbolically in Python using CasADi [57] and exported the resulting OCP to acados [55], which generates and compiles highly optimized C code via its code-generation pipeline. The solver uses the acados SQP-RTI scheme, ERK-4 integration, and the interior-point QP solver HPIPM [56], and a Gauss-Newton Hessian approximation of the OCP. Simulations were performed on the compute nodes of a throughput cluster, each with two 16-core $2.90\,\mathrm{GHz}$ Intel Xeon Gold 6326 (Ice Lake) processors, i.e. 32 physical cores, and $250\,\mathrm{GB}$ of RAM. Every run was pinned to a single core, so the per-robot solve times reported in § 3.6 are single-core numbers.

The physics environment was simulated using a custom 2D engine within a fixed $0.8\,\mathrm{m}\times 0.8\,\mathrm{m}$ workspace (horizon $N=20$, control step $\Delta t=0.1\,\mathrm{s}$, up to $K=5$ neighbors). The NMPC and Reynolds penalty weights were tuned via Bayesian optimization with Optuna [58] to maximize swarm coherence while prioritizing safety and solver stability. The lower speed bound $v_{\min}=0.01$ m/s (Table 3) keeps every robot moving instead of letting it stop. The environmental Gaussian light source rotated along a circular trajectory in the phototaxis experiment and remained stationary in the phototactic orbital experiment. All robots were initialized at uniform random positions in a central sub-region covering 90% of the arena $\mathcal{A}$ with a minimum initial separation of $D_{\rm safe}$ and ran for 9000 control steps ($T_{\rm sim} = 900\,\mathrm{s}$), long enough for the swarm to reach and hold a steady state rather than being measured mid-transient.

Metrics and protocol

All reported quantities are means over $10$ random seeds with $95\%$ confidence intervals. Steady-state quantities are averaged over the last $3000$ of the $9000$ steps (the final $300\,\mathrm{s}$). A robot-robot (or robot-obstacle) contact is a distinct overlap of the elliptical footprints, evaluated at the MPC-commanded positions and counted once per overlap episode. All runs include a lightweight reactive collision-response layer (a per-robot heading bias plus a few positional-separation sweeps) matching the deployed system. A contact is detected as soon as the footprints overlap, before the layer reacts. The layer then pushes the overlapping pair apart and biases their headings away from each other, so overlaps stay shallow and brief. Solver feasibility is the fraction of per-robot solves returning a successful acados status. The six formulations share an identical stack, weights, horizon, arena, seed set, and safe distance $D_{\rm safe}$. Only the inter-robot safety and obstacle-avoidance constraints differ. For the look-ahead formulation, the pair $(\alpha, L)$ was selected by a grid screening followed by Bayesian optimization, and for the anchored Slack Dt-CBF, the decay rate $\alpha$ was selected by a direct two-point comparison, each re-verified over the full ten-seed protocol. Formulations (d) and (e) therefore share $\alpha$, which is what makes the look-ahead ablation below a single-variable comparison. The remaining formulations run at the gains inherited from their source implementations rather than at individually optimized ones; a single shared gain would not be meaningful across formulations that differ in structure, and we report this asymmetry because it bounds how the cross-formulation ranking should be read. The complete set of cost weights, CBF gains, slack weights, and environment constants is listed in Appendix § 6,Table 3.

Capacity and Operating Density

We now characterize how densely the swarm can pack in the fixed arena. This section establishes the geometric capacity and selects the operating point used in the later analysis.

Geometric capacity

The swarm size is limited by the safety buffer, not the physical footprint. The elliptical footprint is small, only $\approx 11\%$ of the arena at $N_b=50$, against $\approx 30\%$ when the buffers are also considered. At $D_{\rm safe}=7\,\mathrm{cm}$, the buffers cap the swarm size at an upper limit of $N_{\max}=137$ non-overlapping safety disks of radius $D_{\rm safe}/2$ (Table 1). This packing results in $\approx 82\%$ area coverage, just below the ideal hexagonal $90.7\%$ because of boundary and finite-size effects. The ceiling depends on the safety distance, ranging from $N_{\rm max}=188$ at $D_{\rm safe}=6\,\mathrm{cm}$ to $N_{\rm max}=68$ at $D_{\rm safe}=10\,\mathrm{cm}$. Across these distances, the buffer occupancy at the limit stays near $82$–$83\%$. We therefore use buffer occupancy as a density measure comparable across safety distances.

Operating point

We select the operating point with a closed-loop DMPC sweep over $D_{\rm safe}$ and movement headroom $\Delta N$. The sweep scores occupancy, collisions, feasibility, solve time, and neighbor-buffer utilization. Sweeping over $D_{\rm safe}\in\{7,8\}\,\mathrm{cm}$ over the four barrier and CBF formulations available at the time, it finds both viable, so we set $D_{\rm safe}=D_{\rm safe, obs}=7\,\mathrm{cm}$, the tighter of the two. The working maximum $N_b=100$ then follows from the geometric headroom criterion $N_{\rm op}=N_{\max}-\Delta N$ with $\Delta N=37$ ($\approx 60\%$ buffer occupancy), not from the sweep score. Beyond $N_{\max}$ the safety buffers can no longer be packed without forced overlap, and the swarm stays feasible only through the slacked barrier.

Safety Formulation Comparison

The inter-robot and obstacle safety constraints can be formulated in several ways. To map this design space, we compare the six formulations (a)–(f) of § 2.5.3 and § 2.5.4 for the phototaxis mode in the absence of obstacles. The slacked variants carry the feasibility-recovery idea of [39] over to our distributed swarm. Under the identical setup of § 3.1 (same stack, solver, and closed-loop simulator), the six formulations differ only in how they enforce inter-robot and obstacle safety, so the comparison isolates that single design choice. Each configuration is evaluated at $N_b \in \{10,25,50,75,100,125,137\}$ over 10 seeds, with $N_{\max}=137$ marked by the dashed vertical line in the figures. Figure 3 reports the solver feasibility, the collisions, and the mean QP solve time versus swarm size $N_b$, and Table 2 summarizes the low-density operating point $N_b = 25$ alongside a high-density operating point $N_b = 100$.

Feasibility

Solver feasibility separates the hard and slacked formulations most sharply (Figure 3a). The two hard formulations hold only about $90\%$ to $92\%$ feasibility at $N_b = 25$. When a solver is infeasible (roughly $8\%$ to $10\%$ of the per-robot OCPs), the robot falls back to a reactive potential-field controller (goal attraction with inter-robot, obstacle, and wall repulsion). This matches the recursive-infeasibility limitation that [38] identify for the hard discrete-time CBF. Softening the barrier with slack variables removes the problem. All four slacked formulations (Slack BF, Slack Dt-CBF, LA Slack Dt-CBF, and Slack Dt-HOCBF) hold at least $99.4\%$ feasibility across the entire density range, and $100\%$ up to $N_b=25$; the minimum over all densities and all four is $99.47\%$, at the packing limit. The small remaining infeasibility comes from the hard arena-boundary constraint, not the safety constraint. At high density, robots pressed against the arena walls can briefly make it infeasible. Feasibility alone, however, does not separate these four. As the collision results show next, a formulation can be feasible everywhere yet still fail to keep robots apart.

Collisions

The formulations diverge in collision counts, not feasibility (Figure 3b), and they separate by structure rather than by barrier order. At $N_b=100$ the three formulations that pair a barrier anchor with a rate condition (d)–(f) average $4.2\times10^{4}$ contacts, against $1.02\times10^{5}$ for the three without an anchor (a)–(c), a factor of $2.4$. The two groups do not overlap at any density in the sweep. Carrying an anchor alongside the rate condition therefore matters more than which rate condition is used.

Two further results follow within the anchored group. First, the look-ahead handle alone reduces contacts by a factor of $2.5$ at low and moderate density, falling to $1.2$ at the packing limit, where it takes them from $1.05\times10^{5}$ to $8.6\times10^{4}$. Formulations (d) and (e) are identical apart from the handle, at the same decay rate $\alpha=0.9$, with disjoint confidence intervals at every density $N_b\ge50$ and no measurable solve-time cost. The reduction is therefore due to the handle, which gives the angular velocity direct authority over the constraint so a robot can steer around a neighbor rather than only brake. Second, the higher-order chain helps at moderate density and stops helping at the packing limit. The Slack Dt-HOCBF is second best in small swarms, but by $N_b=137$ it falls behind the LA Slack Dt-CBF, $1.19\times10^{5}$ against $8.6\times10^{4}$. It does not fail by becoming infeasible, holding $100\%$ feasibility across the whole range against $99.5\%$ for (d) and (e). What it gives up is margin, and not just at the packing limit. The chain lets its robots run closer than the look-ahead variant at every density (Figure 4a), $5.02$ against $5.31\,\mathrm{cm}$ at the packing limit. Its two-step condition allows the barrier to shrink further before the robot has to act. That costs nothing when the arena has unoccupied space. As the available area shrinks, a smaller gap leads to contact. On this task, the deployed LA Slack Dt-CBF is best or tied-best at every density except $N_b=50$, where it is second behind the Slack Dt-HOCBF with disjoint confidence intervals. Elsewhere, the look-ahead advantage over the chain widens as the arena fills, from $1.11\times$ at $N_b=100$ to $1.39\times$ at $N_b=137$.

Safety margins

The formulations also differ in the safety margin they sustain (Figure 4). These runs use the deployed closed loop with the reactive layer active (§ 3.1), and are measured on post-steering positions, so they never show interpenetration. Each value is a mean over robots and steps of the nearest-neighbor distance, not a per-step minimum. The margins separate along the same structural line as the collision counts. On the inter-robot axis, the three anchored formulations (d)–(f) hold $10.1$–$10.6\,\mathrm{cm}$ at $N_b=10$ and $5.0$–$5.3\,\mathrm{cm}$ at the packing limit. The three without an anchor (a)–(c) fall from $7.8$–$8.6\,\mathrm{cm}$ to $4.3$–$4.5\,\mathrm{cm}$. Comparing group means, the anchor sustains roughly $0.7\,\mathrm{cm}$ of extra clearance at maximum crowding. The deployed LA Slack Dt-CBF is widest of all at high density, with $5.82\,\mathrm{cm}$ at $N_b=100$ and $5.31\,\mathrm{cm}$ at $N_{\max}=137$. At $N_b=25$, only the anchored formulations still hold a clear margin above $D_{\rm safe}$, at $7.6$–$8.8\,\mathrm{cm}$. The unanchored group sits at $6.8$–$7.1\,\mathrm{cm}$, at or just either side of it. Beyond that density, every formulation operates inside the nominal buffer, already $6.5$–$6.8$ against $5.8$–$5.9\,\mathrm{cm}$ at $N_b=50$. The task drives this, not the choice of formulation. Phototaxis gathers the swarm at the light peak, so the robots there sit far closer than the arena average of $30\%$ buffer coverage at $N_b=50$ (Table 1) would suggest, and neighboring buffers have to overlap. What the anchor leads to is a consistently larger share of the space that remains, and from $N_b=75$ upward the formulation that holds the widest spacing also records the fewest contacts.

On the robot-obstacle axis, the ordering is similar, but the effect concentrates at the low-density end. At $N_b=10$ the Slack Dt-CBF and the LA Slack Dt-CBF (d),(e) keep about $5.6\,\mathrm{cm}$ of clearance each, against $1.8$–$3.8\,\mathrm{cm}$ for the remaining four, the Slack Dt-HOCBF included. From $N_b=100$ and up, every formulation converges to $1.6$–$1.7\,\mathrm{cm}$. At that density, robots are pressed against obstacles by their neighbors rather than by their own control choices, so the safety formulation can no longer set the obstacle clearance. This is the margin-side counterpart of the collision result that obstacles stop mattering once inter-robot congestion dominates (§ 3.5).

Computational cost

The LA Slack Dt-CBF costs slightly more to solve because of its added look-ahead constraints. Its per-robot solve time is about $3.7\,\mathrm{ms}$ at $N_b = 25$, roughly $1.0\,\mathrm{ms}$ above the cheapest formulations (Figure 3c). Each robot optimizes against a fixed number of neighbors, so this per-robot cost barely changes with swarm size and stays within $3.66$–$3.74\,\mathrm{ms}$ across the sweep. It is therefore well inside the $100\,\mathrm{ms}$ budget of a $10\,\mathrm{Hz}$ loop at every density. § 3.6 reports the distribution and separates this per-robot cost from the per-step wall time of the single-process simulation, which does grow with $N_b$. The extra cost is small, so on phototaxis, the deployment target, the LA Slack Dt-CBF gives the best trade-off across feasibility, collisions, and safety margin. Therefore, we adopt it as the deployed formulation for the collective behavior simulations.

Emergence of Collective Behavior

We characterize the emergent behavior of the swarm under the DMPC framework. The LA Slack Dt-CBF formulation governs both inter-robot separation and obstacle avoidance. We study three tasks, namely coordinated phototaxis, phototactic orbiting, and cohesion without a stimulus. Figure 5 shows representative snapshots at $N_b=50$. The snapshots show how the local coupling terms regulate spatial density and heading agreement while each agent pursues its objective, without global communication or a central coordinator. With two and with twelve static obstacles in the arena, every behavior still holds its form while the swarm moves around the obstacles (Figure 8 and Figure 9 in Appendix § 5).

We assess coordination among the robots using two distinct metrics. The first is the polar order parameter calculated locally over only the set of active neighbors,

where $\mathfrak{N}_{l}$ is the set of active neighbors of each robot $l$ inclusive of the robot $l$, and $K_l = |\mathfrak{N}_l|$ is the number of active neighbors for $l$. When all the robots in every neighborhood are in near-perfect alignment, $\Phi_{\rm loc} \approx 1$, whereas in a disordered state with random robot orientations, $\Phi_{\rm loc} \approx 0$. This metric is widely used in biophysics to quantify order in a moving flock or swarm [8, 3]. The second metric is the nearest-neighbor distance averaged over all robots,

where $k$ is the nearest neighbor of robot $i$. This metric finds pervasive use in quantifying internal group structure in natural and artificial flocks [64, 63]. While the NND varies with animal species, it unequivocally reaches a stable value at long time scales. The top and bottom rows of Figure 6 elucidate the variation in the local polar order and long-time NND for different swarm sizes and different obstacle numbers, under distinct modes of operation.

Phototaxis Behavior

We first examine the full DMPC formulation with all Reynolds steering rules (§ 2.5) active, demonstrating gradient-following behavior, followed by isolating the contribution of the social coupling by removing these terms entirely.

With Reynolds Flocking Rules – Collective Behavior

In phototaxis mode, each robot ascends the light gradient while the cohesion and alignment costs, together with the CBF separation constraints, produce collective motion. The CBF constrains the closing velocity toward each neighbor’s predicted position, so avoidance starts early, without the over-braking or deadlocks of purely reactive methods [53]. The swarm stays compact and heading-aligned as it tracks the rotating source (Figure 5a). The local polar order reduces with increasing swarm size, as radial phototaxis uses a convergent field that brings agents together from all directions, producing head-on encounters that locally cancel heading agreement (Figure 6a). A cluttered environment further worsen the alignment by introducing additional collision avoidance, which de-aligns the local heading. The NND (Figure 6d), on the other hand, decreases as the space for movement in the domain decreases (increasing swarm size), partly also attributable to the nature of the phototaxis mode, which attempts to converge the robots to a small location. Beyond $N_b=50$, the NND falls below the safe buffer distance of $D_{\rm safe} = 7\, \mathrm{cm}$, indicating some amount of buffer overlap. As the swarm size approaches the geometric packing limit, the phototaxis gradient forces the robots into oblique, blocking contacts, causing NND to fall as low as $5\,\mathrm{cm}$.

Without Reynolds Flocking Rules – Swarm Collapse

In the absence of the cohesion and alignment costs and the separation constraints, the collective behavior disappears. Each robot climbs the gradient on its own, so the swarm arrives at the light peak without spacing or heading agreement (Figure 5b). The swarm effectively collapses onto a small area in the domain, clearly attributing the occurrence of flocking to the Reynolds steering costs.

Phototactic Orbital Behavior

In the phototactic orbital mode, each robot follows the tangent of the sensed intensity level set through the gradient-rotation reference (9), orbiting the source without global localization. The CBF constraints (26) prevent collisions as robots compete for the same orbit while the cohesion and alignment costs keep the swarm compact and co-rotating. From arbitrary initial positions, it converges to a stable orbital formation (Figure 5c). Although the phototactic orbital mode is susceptible to crowding similar to the phototaxis mode, the decay of the local polar order is more muted in the former (Figure 6b). This is a direct consequence of the co-location of local neighbors along concentric orbital paths, leading to better heading agreement, even for large swarms. Similarly, the NND also settles slightly below the safe buffer distance even for $N_b=137$, as the locally tangential motion along the orbital causes the robots to align along their major axes, allowing for sufficient separation among them (Figure 6e). Both metrics are only mildly affected by the addition of obstacles to the system.

Flocking in the Absence of Stimuli

In this mode, the light field is disabled, providing the robots no global stimulus. Hence, each robot holds its heading with a small random turn each step. The swarm drifts outward and settles into a band along the walls, with an empty interior (Figure 5d). Once a robot reaches a wall, it stays there until its heading drifts away, so the walls fill faster than they empty. Cohesion and alignment order the band, so the swarm circulates the arena similar to an edge current [62].

The flock remains relatively well-aligned at low swarm sizes, with heading disagreements as the domain starts getting crowded with robots and obstacles due to frequent steering for collision avoidance (Figure 6c). This is an expected outcome, showing that the alignment and cohesion costs work as intended. The band structure resembles the phototactic orbital, despite the difference in shape of the orbit. This is evident from the NND (Figure 6f), which follows a similar trend as in the phototactic orbital mode.

Collision Behavior Across Density and Tasks

The collision load of the deployed LA Slack Dt-CBF depends on both swarm size and task. We measure both here, across three obstacle environments.

Collisions versus swarm size

We sweep $N_b\in\{10,25,50,75,100,125,137\}$ across three environments, namely obstacle-free, few-obstacle ($2$), and cluttered ($12$, positions randomized per seed). We count distinct robot-robot contacts over the $9000$-step run (Figure 7). Contacts stay negligible in small swarms. Free space averages $1.0$ contact at $N_b=25$, with the nearest neighbor $8.8\,\mathrm{cm}$ away on average, well above $D_{\rm safe}$. Contacts then climb steeply, from $4.2\times10^{2}$ at $N_b=50$ to $8.6\times10^{4}$ at the packing limit, where the mean nearest-neighbor separation has fallen to $5.3\,\mathrm{cm}$. The growth is strongly linear. Contacts rise by a factor of $87$ between $N_b=50$ and $N_b=100$ while the robot count only doubles. Near the packing limit the full $D_{\rm safe}$ margin cannot be held for every pair, so the slacked CBF admits transient contacts while the solver stays feasible (§ 3.3).

Obstacles dominate in small swarms and hardly matter in large ones (Figure 7). At $N_b=25$ the cluttered arena produces three orders of magnitude more contacts than free space. By $N_b=100$ that gap has closed to a factor of $1.6$, and by $N_b=125$ to $1.3$. The crossover falls between $N_b=50$ and $N_b=75$. The environment therefore dominates in small swarms, while in large swarms collisions come from the packing density.

Effectiveness of avoidance is set by task geometry

Considering an identical controller, the number of contacts spans two orders of magnitude across the three behaviors, a factor of $103$ between the best and the worst at $N_b=100$ (Figure 7). Phototaxis concentrates the swarm at the peak intensity of the light source and has the highest number of contacts at $3.7\times10^{4}$. The orbital mode spreads robots along an intensity level set and has far fewer contacts at $3.5\times10^{2}$. Cohesion without a stimulus falls in between at $2.7\times10^{3}$. However, this order of effectiveness is not preserved across different swarm sizes. The stimulus-free mode is the worst in small swarm, with $25$ contacts at $N_b=10$, whereas the other modes record none. The cause is how robots meet, not how crowded they are. At $N_b=10$ the stimulus-free swarm is the most spread out of the three, with a mean nearest-neighbor distance of $16.1\,\mathrm{cm}$ against $9.9\,\mathrm{cm}$ for phototaxis. The difference is in their headings. A shared gradient gives phototaxis and orbiting a common direction, so two robots that come close are already traveling the same way. Without a stimulus each robot follows its own random heading, and only the local coupling keeps neighbors together. In large swarms the same absence of a shared target protects it, since the robots never concentrate. At every larger swarm size orbiting stays safest and phototaxis worst, except in the cluttered arena where the stimulus-free mode moves ahead of orbiting at $N_b=50$ and $N_b=75$. Obstacles raise contacts in every mode, but the orbital mode pays the largest penalty, $14.5\times$ at $N_b=100$ against $1.6\times$ for phototaxis. Obstacles sit directly on the orbit and force detours. The safety formulation decides how large a swarm can grow before contacts appear. The task decides how many contacts it has at any given size.