Machine Learning how to Industrial AI & Operational Intelligence Simulating Drones: A Boids Flocking Algorithm in Python

Simulating Drones: A Boids Flocking Algorithm in Python

How simple local rules create realistic drone swarm behavior. In this article, we build a 2D Boids-style flocking simulation step by step in Python.

1. From Birds to Drone Swarms

Drone swarms are an active research topic in multi‑agent systems, where dozens or hundreds of aerial robots must coordinate without a
centralized controller. A surprisingly effective starting point comes from computer graphics: the Boids flocking algorithm.

Boids models each agent (bird, fish, or drone) as a simple particle that follows a handful of geometric rules. When you simulate many agents
together, complex flocking behavior emerges automatically.

Idea: Instead of planning the full swarm trajectory, we give each drone a few local rules: avoid collisions, align with neighbors, and stay with the group. The swarm behavior appears naturally from these interactions.

2. The Boids Flocking Model

The classic Reynolds Boids model defines three key steering behaviors for each agent based on nearby neighbors.

Separation

Steer away from neighbors that are too close to avoid collisions. Each drone computes a repulsive vector from nearby agents.

Alignment

Steer toward the average heading (velocity direction) of nearby drones so the group moves in a similar direction.

Cohesion

Steer toward the average position of neighbors to keep the swarm from dispersing and maintain a coherent cluster.

World Constraints

Optionally add boundaries or target‑seeking terms (e.g., keep drones inside an operational area or guide them toward a waypoint).

Practical tip: Each rule produces a steering vector. In the update step, we compute a weighted sum of these vectors and add it to the agent’s velocity.

3. Designing a Drone Boid in Python

We will implement a 2D Boids simulation using NumPy and matplotlib for visualization, similar to many teaching examples for flocking behavior in Python.

3.1. Environment and Dependencies

Install the required libraries:

pip install numpy matplotlib

3.2. State Representation

Each drone is a point with position and velocity in 2D. We store them in NumPy arrays:

  • positions: shape (N, 2)
  • velocities: shape (N, 2)
import numpy as np

N_DRONES = 40
WIDTH, HEIGHT = 800, 600

positions = np.random.rand(N_DRONES, 2) * np.array([WIDTH, HEIGHT])
angles = np.random.rand(N_DRONES) * 2 * np.pi
speed = 2.0
velocities = np.stack([np.cos(angles), np.sin(angles)], axis=1) * speed
Why vectors? Using arrays for all drones at once lets us leverage NumPy’s vectorized operations for distance computation and steering forces, which is much faster than pure Python loops for large swarms.

4. Implementing the Three Core Rules

For each drone, we look at neighbors within a perception radius and compute three steering vectors: separation, alignment, and cohesion.

4.1. Neighborhood Detection

def get_neighbors(idx, positions, radius):
    diffs = positions - positions[idx]          # (N, 2)
    dist2 = np.sum(diffs**2, axis=1)
    mask = (dist2 > 0) & (dist2 < radius**2)    # exclude self
    return mask, diffs[mask], dist2[mask]

4.2. Separation: Avoid Collisions

def separation(idx, positions, radius, strength=1.0):
    mask, diffs, dist2 = get_neighbors(idx, positions, radius)
    if not np.any(mask):
        return np.zeros(2)

    # Inverse-square repulsion from neighbors
    directions = -diffs / np.expand_dims(np.sqrt(dist2) + 1e-6, axis=1)
    steer = directions.mean(axis=0)
    return steer * strength

4.3. Alignment: Match Heading

def alignment(idx, positions, velocities, radius, strength=0.05):
    mask, _, _ = get_neighbors(idx, positions, radius)
    if not np.any(mask):
        return np.zeros(2)

    avg_vel = velocities[mask].mean(axis=0)
    steer = avg_vel - velocities[idx]
    return steer * strength

4.4. Cohesion: Move Toward Center of Mass

def cohesion(idx, positions, radius, strength=0.01):
    mask, _, _ = get_neighbors(idx, positions, radius)
    if not np.any(mask):
        return np.zeros(2)

    center = positions[mask].mean(axis=0)
    steer = center - positions[idx]
    return steer * strength

4.5. Boundaries and Speed Limits

def apply_bounds(positions, velocities, width, height, margin=50, turn_factor=0.4):
    for i in range(len(positions)):
        x, y = positions[i]
        if x < margin:
            velocities[i, 0] += turn_factor
        elif x > width - margin:
            velocities[i, 0] -= turn_factor

        if y < margin:
            velocities[i, 1] += turn_factor
        elif y > height - margin:
            velocities[i, 1] -= turn_factor


def limit_speed(velocities, max_speed=4.0):
    speeds = np.linalg.norm(velocities, axis=1, keepdims=True)
    mask = speeds > max_speed
    velocities[mask] = velocities[mask] / speeds[mask] * max_speed
Drone interpretation: In a real swarm, these steering vectors would map to small changes in thrust and heading for each UAV, but the geometric logic is the same.

5. Building the Simulation Loop

We now combine all rules into a single update step and animate the result with matplotlib.animation.FuncAnimation.

from matplotlib import pyplot as plt
from matplotlib import animation


PERCEPTION_RADIUS = 60.0
SEP_RADIUS = 25.0   # smaller radius for separation

def update_boids(positions, velocities, dt=1.0):
    new_vel = velocities.copy()

    for i in range(len(positions)):
        sep = separation(i, positions, SEP_RADIUS, strength=0.8)
        ali = alignment(i, positions, velocities, PERCEPTION_RADIUS, strength=0.08)
        coh = cohesion(i, positions, PERCEPTION_RADIUS, strength=0.02)

        steer = sep + ali + coh
        new_vel[i] += steer * dt

    limit_speed(new_vel, max_speed=4.0)
    apply_bounds(positions, new_vel, WIDTH, HEIGHT)

    positions += new_vel * dt
    return positions, new_vel


# Set up plot
fig, ax = plt.subplots(figsize=(8, 6))
scat = ax.scatter(positions[:, 0], positions[:, 1], s=20, c="tabblue")
ax.set_xlim(0, WIDTH)
ax.set_ylim(0, HEIGHT)
ax.set_xticks([])
ax.set_yticks([])
ax.set_title("Drone Swarm Boids Simulation", fontsize=12)


def animate(frame):
    global positions, velocities
    positions, velocities = update_boids(positions, velocities, dt=1.0)
    scat.set_offsets(positions)
    return scat,


anim = animation.FuncAnimation(
    fig, animate, frames=400, interval=30, blit=True
)

plt.show()
Performance note: The above code is intentionally clear rather than fully optimized. For very large swarms, you can use spatial partitioning (grids, k‑d trees) to avoid the O(N²) neighbor search cost.

6. Full Python Example

The following script is a complete, runnable 2D Boids simulation for drone‑like agents. Save it as drone_boids.py and run with python drone_boids.py.

#!/usr/bin/env python3
"""
Simulating Drones: A Boids Flocking Algorithm in Python

This script implements a simple 2D Boids-style flocking model
to simulate a swarm of drone-like agents using NumPy and matplotlib.
"""

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation

# -----------------------
# Simulation Parameters
# -----------------------
N_DRONES = 60
WIDTH, HEIGHT = 800, 600

BASE_SPEED = 2.0
MAX_SPEED = 4.0

PERCEPTION_RADIUS = 60.0
SEP_RADIUS = 25.0

SEP_STRENGTH = 0.8
ALI_STRENGTH = 0.08
COH_STRENGTH = 0.02

BOUND_MARGIN = 50
BOUND_TURN_FACTOR = 0.4

DT = 1.0


# -----------------------
# Initialization
# -----------------------
def init_swarm(n, width, height, base_speed):
    positions = np.random.rand(n, 2) * np.array([width, height])
    angles = np.random.rand(n) * 2 * np.pi
    velocities = np.stack([np.cos(angles), np.sin(angles)], axis=1) * base_speed
    return positions, velocities


positions, velocities = init_swarm(N_DRONES, WIDTH, HEIGHT, BASE_SPEED)


# -----------------------
# Boids Rules
# -----------------------
def get_neighbors(idx, positions, radius):
    diffs = positions - positions[idx]          # (N, 2)
    dist2 = np.sum(diffs**2, axis=1)
    mask = (dist2 > 0) & (dist2 < radius**2)
    return mask, diffs[mask], dist2[mask]


def separation(idx, positions, radius, strength=1.0):
    mask, diffs, dist2 = get_neighbors(idx, positions, radius)
    if not np.any(mask):
        return np.zeros(2)

    directions = -diffs / (np.sqrt(dist2)[:, None] + 1e-6)
    steer = directions.mean(axis=0)
    return steer * strength


def alignment(idx, positions, velocities, radius, strength=0.05):
    mask, _, _ = get_neighbors(idx, positions, radius)
    if not np.any(mask):
        return np.zeros(2)

    avg_vel = velocities[mask].mean(axis=0)
    steer = avg_vel - velocities[idx]
    return steer * strength


def cohesion(idx, positions, radius, strength=0.01):
    mask, _, _ = get_neighbors(idx, positions, radius)
    if not np.any(mask):
        return np.zeros(2)

    center = positions[mask].mean(axis=0)
    steer = center - positions[idx]
    return steer * strength


def apply_bounds(positions, velocities, width, height,
                 margin=50, turn_factor=0.4):
    for i in range(len(positions)):
        x, y = positions[i]
        if x < margin: velocities[i, 0] += turn_factor elif x > width - margin:
            velocities[i, 0] -= turn_factor

        if y < margin: velocities[i, 1] += turn_factor elif y > height - margin:
            velocities[i, 1] -= turn_factor


def limit_speed(velocities, max_speed):
    speeds = np.linalg.norm(velocities, axis=1, keepdims=True)
    mask = speeds > max_speed
    velocities[mask] = velocities[mask] / speeds[mask] * max_speed


# -----------------------
# Update Step
# -----------------------
def update_boids(positions, velocities, dt=1.0):
    new_vel = velocities.copy()

    for i in range(len(positions)):
        sep = separation(i, positions, SEP_RADIUS, strength=SEP_STRENGTH)
        ali = alignment(i, positions, velocities, PERCEPTION_RADIUS,
                        strength=ALI_STRENGTH)
        coh = cohesion(i, positions, PERCEPTION_RADIUS,
                       strength=COH_STRENGTH)

        steer = sep + ali + coh
        new_vel[i] += steer * dt

    limit_speed(new_vel, MAX_SPEED)
    apply_bounds(positions, new_vel, WIDTH, HEIGHT,
                 margin=BOUND_MARGIN, turn_factor=BOUND_TURN_FACTOR)

    positions += new_vel * dt
    return positions, new_vel


# -----------------------
# Visualization
# -----------------------
fig, ax = plt.subplots(figsize=(8, 6))
scat = ax.scatter(positions[:, 0], positions[:, 1],
                  s=25, c="tab:blue", edgecolors="white", linewidths=0.5)

ax.set_xlim(0, WIDTH)
ax.set_ylim(0, HEIGHT)
ax.set_xticks([])
ax.set_yticks([])
ax.set_title("Drone Swarm Boids Simulation", fontsize=13)


def animate(frame):
    global positions, velocities
    positions, velocities = update_boids(positions, velocities, dt=DT)
    scat.set_offsets(positions)
    return scat,


anim = animation.FuncAnimation(
    fig, animate, frames=600, interval=30, blit=True
)

plt.tight_layout()
plt.show()

7. Extending to Drone Swarms

While this implementation runs in 2D and ignores aerodynamics, it captures the core idea behind decentralized drone swarm control: each
vehicle reacts only to nearby neighbors and simple rules, yet the group as a whole shows coordinated, lifelike motion.

See also  Spiking Neural Networks: A Snntorch Simulation Example

Next Steps

  • Add obstacles and obstacle‑avoidance steering.
  • Promote one or more drones to “leaders” that others tend to follow.
  • Lift the simulation to 3D to better approximate real flight.
  • Couple the algorithm with a physics engine or ROS for real UAVs.
  • Use a genetic algorithm to tune rule weights for specific missions.
Machine learning angle: You can treat the rule weights (separation, alignment, cohesion strengths) as parameters and optimize them with evolutionary strategies or reinforcement learning for search‑and‑rescue, coverage, or formation‑keeping tasks.

Leave a Reply

Your email address will not be published. Required fields are marked *

Related Post