Particle Swarm Optimization

Why Particle Swarm Optimization Matters

Imagine tuning two real-valued settings when evaluating the system is possible but differentiating it is difficult. A single search point can forget useful locations after moving away. Particle swarm optimization, or PSO, instead moves a collection of candidate points while remembering each particle's best visit and the best visit found by the entire swarm.

Each particle has a position and a velocity. Its next velocity combines momentum, attraction toward its own remembered success, and attraction toward a shared success. The particle then moves and is evaluated. Its current position may become worse; the best memories are preserved separately. These memories guide exploration without providing a proof that the global optimum has been found.

This lesson implements a specific synchronous global-best PSO for bounded continuous vectors, with inertia weight, separate per-coordinate random factors, velocity caps, and a defined boundary rule. A shifted sphere supplies an analytic reference solution. You will trace motion, memories, random draws, clipping, and objective calls, then run equivalent Python, C++, and Java programs.

Core Concepts and Terminology

Let particle $i$ at iteration $t$ have position $\mathbf x_i^t\in\mathbb R^d$ and velocity $\mathbf v_i^t$. Its personal best $\mathbf p_i^t$ is the first position attaining the smallest objective it has evaluated so far. The swarm best $\mathbf g^t$ is the first position attaining the smallest objective evaluated by any particle. We minimize a deterministic $f$ inside a box $\ell_j\le x_j\le h_j$.

Before velocity and boundary handling, the inertia-weight update for coordinate $j$ is

$$v_{ij}^{\mathrm{raw}}=w v_{ij}^t+c_1r_{1,ij}^t(p_{ij}^t-x_{ij}^t)+c_2r_{2,ij}^t(g_j^t-x_{ij}^t).$$

The inertia weight $w$ retains momentum; $c_1$ controls personal or cognitive attraction; $c_2$ controls social attraction. The random multipliers are sampled separately for every particle, coordinate, and iteration. In an idealized PSO description they are independent uniforms on $[0,1)$; the reproducible teaching programs use the explicitly defined finite grid $\{0,1/1024,\ldots,1023/1024\}$.

The proposed position uses the new velocity, not the previous one:

$$\mathbf x_i^{\mathrm{raw}}=\mathbf x_i^t+\widehat{\mathbf v}_i,$$

where $\widehat{\mathbf v}_i$ is the component-wise capped velocity. Boundary handling may then modify both position and stored velocity. There is no fixed search step size: velocity depends on momentum, remembered points, random factors, and repair.

The worked objective is $f(x,y)=(x-1)^2+(y+1)^2$ on $[-5,5]^2$, with unique minimum $f(1,-1)=0$. The diagram follows particle $2$ through its first two moves. Its current objective worsens from $1$ to $8.5$, while its personal best remains $1$. That distinction is central to PSO.

How Particle Swarm Optimization Works, Step by Step

The inputs provide initial positions and velocities, box bounds, objective center, per-coordinate velocity caps, seed, iteration limit $K$, coefficients $w,c_1,c_2$, and objective-evaluation budget $B$. The evaluator is the shifted sphere

$$f(\mathbf x)=\sum_{j=0}^{d-1}(x_j-c_j)^2.$$

1. Validate and initialize memories. Copy the input arrays. Require budget for all $N$ initial evaluations; otherwise return EVAL_LIMIT with no evaluated solution and no random draws. Evaluate all initial positions, initialize each personal best with an independent position copy, and choose the earliest minimum as the swarm archive. Record iteration $0$.
2. Check complete-iteration limits. If $K$ iterations are complete, stop with ITER_LIMIT. Otherwise require another $N$ evaluations; if unavailable, stop with EVAL_LIMIT before generating any new motion or random draw.
3. Freeze the social guide. Copy the current swarm archive into a guide shared by every particle in this iteration. All calculations read old positions, velocities, and personal memories. Store next states separately. Improvements found early in the sweep affect the guide only in the next iteration.
4. Draw and compute velocity. For each particle, visit coordinates in increasing order. Draw cognitive multiplier first, then social multiplier. Consume both draws even when a coefficient or attraction difference is zero. Compute inertia, cognitive attraction, and social attraction, then add the first two terms and finally the third.
5. Limit velocity. For a nonnegative cap $V_j$, set

$$\widehat v_{ij}=\min\bigl(V_j,\max(-V_j,v_{ij}^{\mathrm{raw}})\bigr).$$

1. Move and handle the box. Form $z_{ij}=x_{ij}^t+\widehat v_{ij}$. If $z_{ij}$ is strictly outside its interval, clip it to the nearest bound and store velocity $0$ for that coordinate. Otherwise store position $z_{ij}$ and velocity $\widehat v_{ij}$. Landing exactly on a bound is not counted as crossing it; outward velocity is reset only if a later proposal goes outside.
2. Evaluate and update memories. Evaluate the final feasible position once, even if it equals the previous position. Always keep the new current state. Replace a personal best only for a strictly smaller objective; similarly replace the swarm archive only for a strictly smaller objective. Equal scores retain the older memory. Never overwrite a good memory just because the current particle moved.
3. Commit and record. After all particles have moved and been evaluated, commit positions, velocities, scores, and personal memories together, increment the iteration, and store a snapshot. The archive may accumulate records during the sweep, but the already-frozen social guide remains unchanged.

The personal and swarm bests are positions with objective values, not just numeric thresholds. Independent copies keep the stored point consistent with its cached score. Unlike differential evolution's target-wise greedy replacement, this PSO accepts current motion that worsens the objective; only its historical memories have monotonic best scores.

flowchart TD
    A["Validate and copy swarm state"] --> B{"Budget covers initialization?"}
    B -->|No| Z["Return no evaluated point"]
    B -->|Yes| C["Evaluate positions and initialize best memories"]
    C --> D{"Iteration limit reached?"}
    D -->|Yes| Y["Return archive and completed swarm"]
    D -->|No| E{"Budget covers N new evaluations?"}
    E -->|No| Y
    E -->|Yes| F["Freeze social guide and start next state"]
    F --> G["Draw per-coordinate factors and compute velocity"]
    G --> H["Cap velocity, move, and repair box crossings"]
    H --> I["Evaluate current point and retain strict memory improvements"]
    I --> J{"More particles?"}
    J -->|Yes| G
    J -->|No| L["Commit iteration and record snapshot"]
    L --> D

Worked Example with a Full Trace

Use $N=3$ particles in two dimensions, bounds $[-5,5]^2$, center $(1,-1)$, caps $(3,3)$, seed $43$, $K=3$, $w=0.5$, $c_1=c_2=1.5$, and $B=100$. Initial state is:

Particle $1$ wins the initial objective tie because it appears first. All initial personal-best coordinates equal their particle's initial position, so the first iteration has zero cognitive attraction.

The seeded generator uses

$$z_{k+1}=48271z_k\bmod2147483647,\qquad z_0=43.$$

Let $R=2147483646$. For each mapped draw, compute $v=z-1$, reject and retry if $v\ge R-(R\bmod1024)$, and otherwise return $q=v\bmod1024$. Use $r=q/1024$. This rejection mapping balances residue classes under a uniform raw-input model; the seeded generator itself is deterministic, not independent mathematical randomness. The first raw state is $2075653$, so the first mapped draw is $(2075653-1)\bmod1024=4$.

For each particle the draw order is cognitive $x$, social $x$, cognitive $y$, social $y$. The next table gives every mapped draw and every raw/capped velocity. Coordinates, velocities, and scores in the trace tables are rounded to eight decimal places for readability; calculations use the unrounded state recorded by the programs.

For iteration $1$, particle $0$, the frozen guide is $(-1,-4)$. Inertia is $(-1.5,1.5)$ and cognitive attraction is zero. Social attraction is $(1.5(760/1024)3,\;1.5(36/1024)(-8))=(3.33984375,-0.421875)$. Thus raw velocity is $(1.83984375,1.078125)$. Its proposed position has second coordinate $5.078125$, so repair sets that coordinate to $5$ and its stored velocity to $0$.

For particle $2$ in the same iteration, raw velocity $(-5.173828125,-8.90234375)$ is capped to $(-3,-3)$, producing $(0,-1)$ with score $1$. At iteration $2$, that particle still has its personal best and frozen social guide at its current point, so both attractions vanish. Inertia alone moves it to $(-1.5,-2.5)$ with score $8.5$; its personal memory stays at $(0,-1)$ with score $1$.

Here are the proposed positions, feasible positions, stored velocities after boundary handling, and memory decisions for every update:

Each accepted personal record stores an independent copy of that row's feasible point. Iteration $3$ also demonstrates nonzero cognitive attraction for particle $2$: its current point is $(-1.5,-2.5)$, while its personal memory is $(0,-1)$ and its social guide is $(0.0946044921875,-1)$. The inertia, cognitive, and social terms are respectively

$$(-0.75,-0.75),\qquad(0.046142578125,0.263671875),$$

and approximately $(1.1936175227165222,0.134033203125)$. The resulting velocity is approximately $(0.4897601008415222,-0.352294921875)$. Its new current score improves relative to its previous current score $8.5$, but remains worse than personal record $1$; that record is therefore preserved.

The completed-iteration summary is:

At iteration $3$, all current particles are worse than the historical swarm archive. The final best is approximately $(0.09460449,-1)$ with objective $0.81974103$, while the best current objective is $2.48337910$. Returning only the current minimum would discard a better evaluated solution. The analytic optimum remains $(1,-1)$ with value $0$.

For exact generator accounting, these are the end states after each particle and cumulative raw draw counts:

No rejection occurs in this sample: each particle update consumes four raw draws, each iteration consumes twelve, and the run consumes $36$. Final status is ITER_LIMIT, with $12$ objective evaluations, $9$ updates, $9$ moved particles, $6$ strict personal records, $3$ strict swarm records, $3$ velocity-coordinate clips, and $1$ box-coordinate crossing. Swarm records can occur more than once within an iteration; the counter excludes initialization.

Implementation in Python, C++, and Java

Each panel is a complete program with the same inputs, random stream, arithmetic order, boundary rule, and printed output. Use Python $3.10+$, C++17, or Java $17+$. For C++, compile with g++ -std=c++17 -O2 -ffp-contract=off main.cpp -o pso and avoid fast-math transformations. Save the Java panel as Main.java.

Inputs require $1\le N\le32$, $1\le d\le16$, a seed in $1,\ldots,2147483646$, $0\le K\le200$, and $0\le B\le10^6$. Coefficients must be finite with $0\le w\le1$ and $0\le c_1,c_2\le2$. Bounds and center coordinates must be finite and have magnitude at most $10^6$, bounds must be ordered, and every initial position must be inside its box. Velocity caps are finite values in $[0,10^6]$; every initial velocity coordinate must have magnitude at most its cap. Equal position bounds and zero velocity caps are allowed. The center may lie outside the box.

Python accepts bounded integer or float numeric data but rejects booleans; seed, iteration limit, and budget must be actual integers. Inputs must be ordinary finite sequences; Java arrays and rows must be non-null. Invalid parameters raise an exception before searching. The finite ranges keep this sphere's squared differences and the velocity arithmetic finite. Underflow and rounding can still erase tiny steps or distinguishable real-valued objectives.

The search uses binary64 arithmetic, with separate subtraction, multiplication, and addition operations. Each attraction computes coefficient times random multiplier, then multiplies the position difference; inertia plus cognitive attraction is added before social attraction. The objective subtracts, squares, and accumulates coordinates left to right. Java's BigDecimal helper is used only to print the exact double rounded to eight decimal places with half-even rounding; the optimizer itself uses double. Display rounding does not change decisions.

Result contains status, has_value, swarm archive and score, final positions and velocities, current scores, personal-best positions and scores, completed iterations, evaluation/update/move/record/clip/hit counts, generator state/draw/retry counts, snapshots, and per-particle events. A moved particle means at least one position value changed numerically. Clip and boundary counters count coordinates; a strict personal record and a swarm record count particle evaluations, not iterations.

When $B\lt N$, the returned arrays and histories are empty and has_value=false; the zero objective field is a dummy value. With enough initialization budget, a zero-iteration request still evaluates the initial swarm. After initialization, ITER_LIMIT is checked before EVAL_LIMIT. Budget $5$ evaluates the initial three positions but cannot fund another three evaluations; it returns the initial archive without consuming randomness. Repeated and unchanged positions are evaluated again; the program does not silently reuse current evaluations.

The shifted sphere is a teaching evaluator with an independently known minimizer. To use a general deterministic objective, replace sphere and check every returned score is finite before updating memories. Constraints, failed simulations, noise, and changing objectives need explicit policies. The collapsed sample uses three identical positions $(4,4)$ and zero velocities: personal and social differences stay zero forever, preserving objective $34$ even though a feasible point with objective $0$ exists.

from dataclasses import dataclass, field
import math


MOD = 2147483647
RANGE = MOD - 1


class RandomSource:
    def __init__(self, seed):
        self.state, self.draws, self.rejections = seed, 0, 0

    def below(self, bound):
        limit = RANGE - RANGE % bound
        while True:
            self.state = 48271 * self.state % MOD
            self.draws += 1
            value = self.state - 1
            if value < limit:
                return value % bound
            self.rejections += 1


def sphere(point, center):
    total = 0.0
    for x, c in zip(point, center):
        delta = x - c
        term = delta * delta
        total = total + term
    return total


def copy_rows(rows):
    return [row[:] for row in rows]


@dataclass
class Event:
    iteration: int
    particle: int
    cognitive_draws: list
    social_draws: list
    raw_velocity: list
    limited_velocity: list
    raw_position: list
    position: list
    velocity: list
    fitness: float
    personal_improved: bool
    rng_state: int
    rng_draws: int


@dataclass
class Snapshot:
    iteration: int
    positions: list
    velocities: list
    fitness: list
    personal: list
    personal_fitness: list
    best: list
    best_fitness: float
    evaluations: int
    rng_draws: int


@dataclass
class Result:
    status: str = "EVAL_LIMIT"
    has_value: bool = False
    best: list = field(default_factory=list)
    best_fitness: float = 0.0
    positions: list = field(default_factory=list)
    velocities: list = field(default_factory=list)
    fitness: list = field(default_factory=list)
    personal: list = field(default_factory=list)
    personal_fitness: list = field(default_factory=list)
    iterations: int = 0
    evaluations: int = 0
    updates: int = 0
    moved: int = 0
    personal_improvements: int = 0
    global_improvements: int = 0
    velocity_clips: int = 0
    boundary_hits: int = 0
    rng_state: int = 0
    rng_draws: int = 0
    rng_rejections: int = 0
    history: list = field(default_factory=list)
    events: list = field(default_factory=list)


def particle_swarm(initial, initial_velocity, lower, upper, center, velocity_cap,
                   seed, iterations, inertia, cognitive, social, evaluation_budget):
    if any(type(x) is not int for x in [seed, iterations, evaluation_budget]):
        raise ValueError("integer parameters required")
    if not (1 <= seed < MOD and 0 <= iterations <= 200 and
            0 <= evaluation_budget <= 1000000):
        raise ValueError("parameter out of range")

    def number(x):
        return type(x) is int or (type(x) is float and math.isfinite(x))

    if not all(number(x) for x in [inertia, cognitive, social]) or not (
            0 <= inertia <= 1 and 0 <= cognitive <= 2 and 0 <= social <= 2):
        raise ValueError("invalid coefficient")
    dimension, size = len(lower), len(initial)
    if not 1 <= dimension <= 16 or not 1 <= size <= 32:
        raise ValueError("invalid dimension or swarm size")
    if any(len(row) != dimension for row in [upper, center, velocity_cap]):
        raise ValueError("vector size mismatch")
    if any(not number(x) or abs(x) > 1000000 for row in [lower, upper, center] for x in row):
        raise ValueError("invalid bound or center")
    if any(lower[j] > upper[j] for j in range(dimension)):
        raise ValueError("reversed bounds")
    if any(not number(x) or not 0 <= x <= 1000000 for x in velocity_cap):
        raise ValueError("invalid velocity cap")
    if len(initial_velocity) != size:
        raise ValueError("velocity count mismatch")
    for point, velocity in zip(initial, initial_velocity):
        if len(point) != dimension or len(velocity) != dimension:
            raise ValueError("row size mismatch")
        for j in range(dimension):
            if not number(point[j]) or not lower[j] <= point[j] <= upper[j]:
                raise ValueError("invalid initial point")
            if not number(velocity[j]) or abs(velocity[j]) > velocity_cap[j]:
                raise ValueError("invalid initial velocity")
    positions = [list(map(float, row)) for row in initial]
    velocities = [list(map(float, row)) for row in initial_velocity]
    lo, hi, origin, cap = [list(map(float, row)) for row in [lower, upper, center, velocity_cap]]
    inertia, cognitive, social = map(float, [inertia, cognitive, social])
    rng = RandomSource(seed)
    result = Result(rng_state=seed)
    if evaluation_budget < size:
        return result
    fitness = [sphere(row, origin) for row in positions]
    personal, personal_fitness = copy_rows(positions), fitness[:]
    best = min(range(size), key=lambda i: fitness[i])
    result.has_value, result.evaluations = True, size
    result.best, result.best_fitness = positions[best][:], fitness[best]

    def record():
        result.history.append(Snapshot(result.iterations, copy_rows(positions),
            copy_rows(velocities), fitness[:], copy_rows(personal), personal_fitness[:],
            result.best[:], result.best_fitness, result.evaluations, rng.draws))

    record()
    while True:
        if result.iterations == iterations:
            result.status = "ITER_LIMIT"
            break
        if evaluation_budget - result.evaluations < size:
            result.status = "EVAL_LIMIT"
            break
        social_best = result.best[:]
        next_positions, next_velocities, next_fitness = [], [], []
        next_personal, next_personal_fitness = copy_rows(personal), personal_fitness[:]
        for i in range(size):
            draws1, draws2, raw_v, limited_v, raw_x, point, velocity = [], [], [], [], [], [], []
            for j in range(dimension):
                q1, q2 = rng.below(1024), rng.below(1024)
                draws1.append(q1)
                draws2.append(q2)
                r1, r2 = q1 / 1024.0, q2 / 1024.0
                momentum = inertia * velocities[i][j]
                personal_gain = cognitive * r1
                personal_delta = personal[i][j] - positions[i][j]
                attraction1 = personal_gain * personal_delta
                social_gain = social * r2
                social_delta = social_best[j] - positions[i][j]
                attraction2 = social_gain * social_delta
                combined = momentum + attraction1
                value = combined + attraction2
                limited = -cap[j] if value < -cap[j] else cap[j] if value > cap[j] else value
                result.velocity_clips += int(value < -cap[j] or value > cap[j])
                proposed = positions[i][j] + limited
                hit = proposed < lo[j] or proposed > hi[j]
                bounded = lo[j] if proposed < lo[j] else hi[j] if proposed > hi[j] else proposed
                result.boundary_hits += int(hit)
                raw_v.append(value)
                limited_v.append(limited)
                raw_x.append(proposed)
                point.append(bounded)
                velocity.append(0.0 if hit else limited)
            value = sphere(point, origin)
            result.evaluations += 1
            result.updates += 1
            result.moved += int(point != positions[i])
            improved = value < personal_fitness[i]
            if improved:
                next_personal[i], next_personal_fitness[i] = point[:], value
                result.personal_improvements += 1
            if value < result.best_fitness:
                result.best, result.best_fitness = point[:], value
                result.global_improvements += 1
            next_positions.append(point[:])
            next_velocities.append(velocity[:])
            next_fitness.append(value)
            result.events.append(Event(result.iterations + 1, i, draws1, draws2,
                raw_v, limited_v, raw_x, point, velocity, value, improved, rng.state, rng.draws))
        positions, velocities, fitness = next_positions, next_velocities, next_fitness
        personal, personal_fitness = next_personal, next_personal_fitness
        result.iterations += 1
        record()
    result.positions, result.velocities = copy_rows(positions), copy_rows(velocities)
    result.fitness = fitness[:]
    result.personal, result.personal_fitness = copy_rows(personal), personal_fitness[:]
    result.rng_state, result.rng_draws = rng.state, rng.draws
    result.rng_rejections = rng.rejections
    return result


def vector(point):
    return "(" + ",".join(f"{x:.8f}" for x in point) + ")"


def show(label, result):
    print(label)
    for row in result.history:
        print(f"t={row.iteration} best={vector(row.best)} f={row.best_fitness:.8f} "
              f"current_min={min(row.fitness):.8f} eval={row.evaluations} draws={row.rng_draws}")
    best = f"{vector(result.best)}:{result.best_fitness:.8f}" if result.has_value else "NONE"
    print(f"status={result.status} best={best} iterations={result.iterations} "
          f"eval={result.evaluations} updates={result.updates} moved={result.moved} "
          f"personal={result.personal_improvements} global={result.global_improvements} "
          f"vclips={result.velocity_clips} hits={result.boundary_hits} "
          f"draws={result.rng_draws} rejections={result.rng_rejections} rng={result.rng_state}")


def main():
    points, velocities = [[-4, 4], [-1, -4], [3, 2]], [[-3, 3], [0.5, 0], [0, -1]]
    for label, budget in [("sphere", 100), ("budget", 5), ("zero-budget", 0)]:
        show(label, particle_swarm(points, velocities, [-5, -5], [5, 5], [1, -1],
            [3, 3], 43, 3, .5, 1.5, 1.5, budget))
    show("collapsed", particle_swarm([[4, 4]] * 3, [[0, 0]] * 3, [-5, -5],
        [5, 5], [1, -1], [3, 3], 7, 2, .5, 1.5, 1.5, 100))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
using namespace std;
using I64 = int64_t;
using Rows = vector<vector<double>>;
constexpr I64 MOD = 2147483647, RANGE = MOD - 1;
struct RandomSource {
    I64 state; int draws = 0, rejections = 0;
    explicit RandomSource(I64 seed) : state(seed) {}
    int below(int bound) {
        I64 limit = RANGE - RANGE % bound;
        while (true) {
            state = 48271 * state % MOD; ++draws;
            I64 value = state - 1;
            if (value < limit) return int(value % bound);
            ++rejections;
        }
    }
};
double sphere(const vector<double>& point, const vector<double>& center) {
    double total = 0.0;
    for (size_t j = 0; j < point.size(); ++j) {
        double delta = point[j] - center[j]; double term = delta * delta;
        total = total + term;
    }
    return total;
}
struct Event {
    int iteration, particle;
    vector<int> cognitive_draws, social_draws;
    vector<double> raw_velocity, limited_velocity, raw_position, position, velocity;
    double fitness; bool personal_improved; I64 rng_state; int rng_draws;
};
struct Snapshot {
    int iteration; Rows positions, velocities; vector<double> fitness;
    Rows personal; vector<double> personal_fitness, best; double best_fitness;
    int evaluations, rng_draws;
};
struct Result {
    string status = "EVAL_LIMIT"; bool has_value = false;
    vector<double> best; double best_fitness = 0.0;
    Rows positions, velocities; vector<double> fitness;
    Rows personal; vector<double> personal_fitness;
    int iterations = 0, evaluations = 0, updates = 0, moved = 0;
    int personal_improvements = 0, global_improvements = 0, velocity_clips = 0, boundary_hits = 0;
    I64 rng_state = 0; int rng_draws = 0, rng_rejections = 0;
    vector<Snapshot> history; vector<Event> events;
};
Result particle_swarm(const Rows& initial, const Rows& initial_velocity,
        const vector<double>& lower, const vector<double>& upper, const vector<double>& center,
        const vector<double>& velocity_cap, I64 seed, int iterations, double inertia,
        double cognitive, double social, int evaluation_budget) {
    if (seed < 1 || seed >= MOD || iterations < 0 || iterations > 200 || evaluation_budget < 0 || evaluation_budget > 1000000)
        throw invalid_argument("parameter out of range");
    if (!isfinite(inertia) || !isfinite(cognitive) || !isfinite(social) || inertia < 0 || inertia > 1 ||
        cognitive < 0 || cognitive > 2 || social < 0 || social > 2)
        throw invalid_argument("invalid coefficient");
    int dimension = int(lower.size()), size = int(initial.size());
    if (dimension < 1 || dimension > 16 || size < 1 || size > 32)
        throw invalid_argument("invalid dimension or swarm size");
    if (upper.size() != lower.size() || center.size() != lower.size() || velocity_cap.size() != lower.size())
        throw invalid_argument("vector size mismatch");
    for (int j = 0; j < dimension; ++j) {
        if (!isfinite(lower[j]) || !isfinite(upper[j]) || !isfinite(center[j]) ||
            abs(lower[j]) > 1000000 || abs(upper[j]) > 1000000 || abs(center[j]) > 1000000)
            throw invalid_argument("invalid bound or center");
        if (lower[j] > upper[j]) throw invalid_argument("reversed bounds");
        if (!isfinite(velocity_cap[j]) || velocity_cap[j] < 0 || velocity_cap[j] > 1000000)
            throw invalid_argument("invalid velocity cap");
    }
    if (initial_velocity.size() != initial.size()) throw invalid_argument("velocity count mismatch");
    for (int i = 0; i < size; ++i) {
        if (initial[i].size() != lower.size() || initial_velocity[i].size() != lower.size())
            throw invalid_argument("row size mismatch");
        for (int j = 0; j < dimension; ++j) {
            if (!isfinite(initial[i][j]) || initial[i][j] < lower[j] || initial[i][j] > upper[j])
                throw invalid_argument("invalid initial point");
            if (!isfinite(initial_velocity[i][j]) || abs(initial_velocity[i][j]) > velocity_cap[j])
                throw invalid_argument("invalid initial velocity");
        }
    }
    auto positions = initial, velocities = initial_velocity;
    auto lo = lower, hi = upper, origin = center, cap = velocity_cap;
    RandomSource rng(seed); Result result; result.rng_state = seed;
    if (evaluation_budget < size) return result;
    vector<double> fitness;
    for (const auto& row : positions) fitness.push_back(sphere(row, origin));
    auto personal = positions; auto personal_fitness = fitness;
    int best = int(min_element(fitness.begin(), fitness.end()) - fitness.begin());
    result.has_value = true; result.evaluations = size;
    result.best = positions[best]; result.best_fitness = fitness[best];
    auto record = [&]() {
        result.history.push_back(Snapshot{result.iterations, positions, velocities, fitness,
            personal, personal_fitness, result.best, result.best_fitness, result.evaluations, rng.draws});
    };
    record();
    while (true) {
        if (result.iterations == iterations) { result.status = "ITER_LIMIT"; break; }
        if (evaluation_budget - result.evaluations < size) { result.status = "EVAL_LIMIT"; break; }
        auto social_best = result.best;
        Rows next_positions, next_velocities; vector<double> next_fitness;
        auto next_personal = personal; auto next_personal_fitness = personal_fitness;
        for (int i = 0; i < size; ++i) {
            vector<int> draws1, draws2;
            vector<double> raw_v, limited_v, raw_x, point, velocity;
            for (int j = 0; j < dimension; ++j) {
                int q1 = rng.below(1024), q2 = rng.below(1024);
                draws1.push_back(q1); draws2.push_back(q2);
                double r1 = q1 / 1024.0, r2 = q2 / 1024.0;
                double momentum = inertia * velocities[i][j];
                double personal_gain = cognitive * r1;
                double personal_delta = personal[i][j] - positions[i][j];
                double attraction1 = personal_gain * personal_delta;
                double social_gain = social * r2;
                double social_delta = social_best[j] - positions[i][j];
                double attraction2 = social_gain * social_delta;
                double combined = momentum + attraction1;
                double value = combined + attraction2;
                double limited = value < -cap[j] ? -cap[j] : value > cap[j] ? cap[j] : value;
                result.velocity_clips += value < -cap[j] || value > cap[j];
                double proposed = positions[i][j] + limited;
                bool hit = proposed < lo[j] || proposed > hi[j];
                double bounded = proposed < lo[j] ? lo[j] : proposed > hi[j] ? hi[j] : proposed;
                result.boundary_hits += hit;
                raw_v.push_back(value); limited_v.push_back(limited); raw_x.push_back(proposed);
                point.push_back(bounded); velocity.push_back(hit ? 0.0 : limited);
            }
            double value = sphere(point, origin); ++result.evaluations; ++result.updates;
            result.moved += point != positions[i];
            bool improved = value < personal_fitness[i];
            if (improved) {
                next_personal[i] = point; next_personal_fitness[i] = value; ++result.personal_improvements;
            }
            if (value < result.best_fitness) {
                result.best = point; result.best_fitness = value; ++result.global_improvements;
            }
            next_positions.push_back(point); next_velocities.push_back(velocity); next_fitness.push_back(value);
            result.events.push_back(Event{result.iterations + 1, i, draws1, draws2, raw_v,
                limited_v, raw_x, point, velocity, value, improved, rng.state, rng.draws});
        }
        positions = move(next_positions); velocities = move(next_velocities); fitness = move(next_fitness);
        personal = move(next_personal); personal_fitness = move(next_personal_fitness);
        ++result.iterations; record();
    }
    result.positions = positions; result.velocities = velocities; result.fitness = fitness;
    result.personal = personal; result.personal_fitness = personal_fitness;
    result.rng_state = rng.state; result.rng_draws = rng.draws; result.rng_rejections = rng.rejections;
    return result;
}
string vector_text(const vector<double>& point) {
    ostringstream out; out << fixed << setprecision(8) << '(';
    for (size_t j = 0; j < point.size(); ++j) { if (j) out << ','; out << point[j]; }
    out << ')'; return out.str();
}
void show(const string& label, const Result& result) {
    cout << fixed << setprecision(8) << label << '\n';
    for (const auto& row : result.history)
        cout << "t=" << row.iteration << " best=" << vector_text(row.best) << " f=" << row.best_fitness
             << " current_min=" << *min_element(row.fitness.begin(), row.fitness.end())
             << " eval=" << row.evaluations << " draws=" << row.rng_draws << '\n';
    cout << "status=" << result.status << " best=";
    if (result.has_value) cout << vector_text(result.best) << ':' << result.best_fitness;
    else cout << "NONE";
    cout << " iterations=" << result.iterations << " eval=" << result.evaluations << " updates=" << result.updates
         << " moved=" << result.moved << " personal=" << result.personal_improvements << " global=" << result.global_improvements
         << " vclips=" << result.velocity_clips << " hits=" << result.boundary_hits << " draws=" << result.rng_draws
         << " rejections=" << result.rng_rejections << " rng=" << result.rng_state << '\n';
}
int main() {
    Rows points, velocities;
    points.push_back(vector<double>{-4, 4}); points.push_back(vector<double>{-1, -4}); points.push_back(vector<double>{3, 2});
    velocities.push_back(vector<double>{-3, 3}); velocities.push_back(vector<double>{.5, 0}); velocities.push_back(vector<double>{0, -1});
    vector<double> lo{-5, -5}, hi{5, 5}, center{1, -1}, cap{3, 3};
    show("sphere", particle_swarm(points, velocities, lo, hi, center, cap, 43, 3, .5, 1.5, 1.5, 100));
    show("budget", particle_swarm(points, velocities, lo, hi, center, cap, 43, 3, .5, 1.5, 1.5, 5));
    show("zero-budget", particle_swarm(points, velocities, lo, hi, center, cap, 43, 3, .5, 1.5, 1.5, 0));
    show("collapsed", particle_swarm(Rows(3, vector<double>{4, 4}), Rows(3, vector<double>{0, 0}),
         lo, hi, center, cap, 7, 2, .5, 1.5, 1.5, 100));
    return 0;
}

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;

public class Main {
    static final long MOD=2147483647L,RANGE=MOD-1;
    static class RandomSource {
        long state;int draws,rejections;
        RandomSource(long seed) { state=seed; }
        int below(int bound) {
            long limit=RANGE-RANGE%bound;
            while (true) {
                state=48271*state%MOD;++draws;long value=state-1;
                if (value<limit) return (int)(value%bound);
                ++rejections;
            }
        }
    }
    static double sphere(double[] point,double[] center) {
        double total=0.0;
        for (int j=0;j<point.length;++j) {
            double delta=point[j]-center[j];double term=delta*delta;total=total+term;
        }
        return total;
    }
    static double[][] copy(double[][] rows) {
        double[][] out=new double[rows.length][];
        for (int i=0;i<rows.length;++i) out[i]=rows[i].clone();return out;
    }
    static class Event {
        int iteration,particle;int[] cognitive_draws,social_draws;
        double[] raw_velocity,limited_velocity,raw_position,position,velocity;
        double fitness;boolean personal_improved;long rng_state;int rng_draws;
        Event(int t,int i,int[] q1,int[] q2,double[] rv,double[] lv,double[] rx,
              double[] x,double[] v,double f,boolean ok,long rs,int rd) {
            iteration=t;particle=i;cognitive_draws=q1;social_draws=q2;raw_velocity=rv;
            limited_velocity=lv;raw_position=rx;position=x;velocity=v;fitness=f;
            personal_improved=ok;rng_state=rs;rng_draws=rd;
        }
    }
    static class Snapshot {
        int iteration;double[][] positions,velocities;double[] fitness;
        double[][] personal;double[] personal_fitness,best;double best_fitness;int evaluations,rng_draws;
        Snapshot(int t,double[][] x,double[][] v,double[] f,double[][] p,double[] pf,double[] b,double bf,int ev,int rd) {
            iteration=t;positions=copy(x);velocities=copy(v);fitness=f.clone();personal=copy(p);
            personal_fitness=pf.clone();best=b.clone();best_fitness=bf;evaluations=ev;rng_draws=rd;
        }
    }
    static class Result {
        String status="EVAL_LIMIT";boolean has_value;double[] best=new double[0];double best_fitness;
        double[][] positions=new double[0][],velocities=new double[0][];double[] fitness=new double[0];
        double[][] personal=new double[0][];double[] personal_fitness=new double[0];
        int iterations,evaluations,updates,moved,personal_improvements,global_improvements,velocity_clips,boundary_hits;
        long rng_state;int rng_draws,rng_rejections;
        List<Snapshot> history=new ArrayList<>();List<Event> events=new ArrayList<>();
    }
    static Result particle_swarm(double[][] initial,double[][] initial_velocity,double[] lower,
            double[] upper,double[] center,double[] velocity_cap,long seed,int iterations,
            double inertia,double cognitive,double social,int evaluation_budget) {
        if (seed<1 || seed>=MOD || iterations<0 || iterations>200 || evaluation_budget<0 || evaluation_budget>1000000)
            throw new IllegalArgumentException("parameter out of range");
        if (!Double.isFinite(inertia) || !Double.isFinite(cognitive) || !Double.isFinite(social) ||
            inertia<0 || inertia>1 || cognitive<0 || cognitive>2 || social<0 || social>2)
            throw new IllegalArgumentException("invalid coefficient");
        int dimension=lower.length,size=initial.length;
        if (dimension<1 || dimension>16 || size<1 || size>32)
            throw new IllegalArgumentException("invalid dimension or swarm size");
        if (upper.length!=dimension || center.length!=dimension || velocity_cap.length!=dimension)
            throw new IllegalArgumentException("vector size mismatch");
        for (int j=0;j<dimension;++j) {
            if (!Double.isFinite(lower[j]) || !Double.isFinite(upper[j]) || !Double.isFinite(center[j]) ||
                Math.abs(lower[j])>1000000 || Math.abs(upper[j])>1000000 || Math.abs(center[j])>1000000)
                throw new IllegalArgumentException("invalid bound or center");
            if (lower[j]>upper[j]) throw new IllegalArgumentException("reversed bounds");
            if (!Double.isFinite(velocity_cap[j]) || velocity_cap[j]<0 || velocity_cap[j]>1000000)
                throw new IllegalArgumentException("invalid velocity cap");
        }
        if (initial_velocity.length!=size) throw new IllegalArgumentException("velocity count mismatch");
        for (int i=0;i<size;++i) {
            if (initial[i].length!=dimension || initial_velocity[i].length!=dimension)
                throw new IllegalArgumentException("row size mismatch");
            for (int j=0;j<dimension;++j) {
                if (!Double.isFinite(initial[i][j]) || initial[i][j]<lower[j] || initial[i][j]>upper[j])
                    throw new IllegalArgumentException("invalid initial point");
                if (!Double.isFinite(initial_velocity[i][j]) || Math.abs(initial_velocity[i][j])>velocity_cap[j])
                    throw new IllegalArgumentException("invalid initial velocity");
            }
        }
        double[][] positions=copy(initial),velocities=copy(initial_velocity);
        double[] lo=lower.clone(),hi=upper.clone(),origin=center.clone(),cap=velocity_cap.clone();
        RandomSource rng=new RandomSource(seed);Result result=new Result();result.rng_state=seed;
        if (evaluation_budget<size) return result;
        double[] fitness=new double[size];int best=0;
        for (int i=0;i<size;++i) {
            fitness[i]=sphere(positions[i],origin);if (fitness[i]<fitness[best]) best=i;
        }
        double[][] personal=copy(positions);double[] personal_fitness=fitness.clone();
        result.has_value=true;result.evaluations=size;result.best=positions[best].clone();result.best_fitness=fitness[best];
        result.history.add(new Snapshot(0,positions,velocities,fitness,personal,personal_fitness,result.best,
            result.best_fitness,result.evaluations,rng.draws));
        while (true) {
            if (result.iterations==iterations) { result.status="ITER_LIMIT";break; }
            if (evaluation_budget-result.evaluations<size) { result.status="EVAL_LIMIT";break; }
            double[] social_best=result.best.clone();
            double[][] next_positions=new double[size][],next_velocities=new double[size][];
            double[] next_fitness=new double[size];
            double[][] next_personal=copy(personal);double[] next_personal_fitness=personal_fitness.clone();
            for (int i=0;i<size;++i) {
                int[] draws1=new int[dimension],draws2=new int[dimension];
                double[] raw_v=new double[dimension],limited_v=new double[dimension],raw_x=new double[dimension];
                double[] point=new double[dimension],velocity=new double[dimension];
                for (int j=0;j<dimension;++j) {
                    int q1=rng.below(1024),q2=rng.below(1024);draws1[j]=q1;draws2[j]=q2;
                    double r1=q1/1024.0,r2=q2/1024.0;
                    double momentum=inertia*velocities[i][j];
                    double personal_gain=cognitive*r1,personal_delta=personal[i][j]-positions[i][j];
                    double attraction1=personal_gain*personal_delta;
                    double social_gain=social*r2,social_delta=social_best[j]-positions[i][j];
                    double attraction2=social_gain*social_delta,combined=momentum+attraction1;
                    double value=combined+attraction2;
                    double limited=value < -cap[j] ? -cap[j] : value>cap[j] ? cap[j] : value;
                    if (value < -cap[j] || value>cap[j]) ++result.velocity_clips;
                    double proposed=positions[i][j]+limited;
                    boolean hit=proposed<lo[j] || proposed>hi[j];
                    double bounded=proposed<lo[j] ? lo[j] : proposed>hi[j] ? hi[j] : proposed;
                    if (hit) ++result.boundary_hits;
                    raw_v[j]=value;limited_v[j]=limited;raw_x[j]=proposed;
                    point[j]=bounded;velocity[j]=hit ? 0.0 : limited;
                }
                double value=sphere(point,origin);++result.evaluations;++result.updates;
                boolean moved=false;for (int j=0;j<dimension;++j) if (point[j]!=positions[i][j]) moved=true;
                if (moved) ++result.moved;
                boolean improved=value<personal_fitness[i];
                if (improved) {
                    next_personal[i]=point.clone();next_personal_fitness[i]=value;++result.personal_improvements;
                }
                if (value<result.best_fitness) {
                    result.best=point.clone();result.best_fitness=value;++result.global_improvements;
                }
                next_positions[i]=point.clone();next_velocities[i]=velocity.clone();next_fitness[i]=value;
                result.events.add(new Event(result.iterations+1,i,draws1,draws2,raw_v,limited_v,raw_x,
                    point,velocity,value,improved,rng.state,rng.draws));
            }
            positions=next_positions;velocities=next_velocities;fitness=next_fitness;
            personal=next_personal;personal_fitness=next_personal_fitness;++result.iterations;
            result.history.add(new Snapshot(result.iterations,positions,velocities,fitness,personal,personal_fitness,
                result.best,result.best_fitness,result.evaluations,rng.draws));
        }
        result.positions=copy(positions);result.velocities=copy(velocities);result.fitness=fitness.clone();
        result.personal=copy(personal);result.personal_fitness=personal_fitness.clone();
        result.rng_state=rng.state;result.rng_draws=rng.draws;result.rng_rejections=rng.rejections;return result;
    }
    static String fixed8(double x) {
        String s=new BigDecimal(x).setScale(8,RoundingMode.HALF_EVEN).toPlainString();
        return Double.doubleToRawLongBits(x)<0 && s.charAt(0)!='-' ? "-"+s : s;
    }
    static String vector(double[] point) {
        StringBuilder s=new StringBuilder("(");
        for (int j=0;j<point.length;++j) { if (j>0) s.append(',');s.append(fixed8(point[j])); }
        return s.append(')').toString();
    }
    static void show(String label,Result result) {
        System.out.println(label);
        for (Snapshot row : result.history) {
            double current=row.fitness[0];for (double x : row.fitness) if (x<current) current=x;
            System.out.println("t="+row.iteration+" best="+vector(row.best)+" f="+fixed8(row.best_fitness)+
                " current_min="+fixed8(current)+" eval="+row.evaluations+" draws="+row.rng_draws);
        }
        String best=result.has_value ? vector(result.best)+":"+fixed8(result.best_fitness) : "NONE";
        System.out.println("status="+result.status+" best="+best+" iterations="+result.iterations+
            " eval="+result.evaluations+" updates="+result.updates+" moved="+result.moved+
            " personal="+result.personal_improvements+" global="+result.global_improvements+
            " vclips="+result.velocity_clips+" hits="+result.boundary_hits+" draws="+result.rng_draws+
            " rejections="+result.rng_rejections+" rng="+result.rng_state);
    }
    public static void main(String[] args) {
        double[][] points=new double[3][],velocities=new double[3][];
        points[0]=new double[]{-4,4};points[1]=new double[]{-1,-4};points[2]=new double[]{3,2};
        velocities[0]=new double[]{-3,3};velocities[1]=new double[]{.5,0};velocities[2]=new double[]{0,-1};
        double[] lo={-5,-5},hi={5,5},center={1,-1},cap={3,3};
        show("sphere",particle_swarm(points,velocities,lo,hi,center,cap,43,3,.5,1.5,1.5,100));
        show("budget",particle_swarm(points,velocities,lo,hi,center,cap,43,3,.5,1.5,1.5,5));
        show("zero-budget",particle_swarm(points,velocities,lo,hi,center,cap,43,3,.5,1.5,1.5,0));
        for (int i=0;i<3;++i) { points[i]=new double[]{4,4};velocities[i]=new double[]{0,0}; }
        show("collapsed",particle_swarm(points,velocities,lo,hi,center,cap,7,2,.5,1.5,1.5,100));
    }
}

Expected standard output, identical in all three languages:

• sphere
• t=0 best=(-1.00000000,-4.00000000) f=13.00000000 current_min=13.00000000 eval=3 draws=0
• t=1 best=(0.00000000,-1.00000000) f=1.00000000 current_min=1.00000000 eval=6 draws=12
• t=2 best=(0.09460449,-1.00000000) f=0.81974103 current_min=0.81974103 eval=9 draws=24
• t=3 best=(0.09460449,-1.00000000) f=0.81974103 current_min=2.48337910 eval=12 draws=36
• status=ITER_LIMIT best=(0.09460449,-1.00000000):0.81974103 iterations=3 eval=12 updates=9 moved=9 personal=6 global=3 vclips=3 hits=1 draws=36 rejections=0 rng=2063162430
• budget
• t=0 best=(-1.00000000,-4.00000000) f=13.00000000 current_min=13.00000000 eval=3 draws=0
• status=EVAL_LIMIT best=(-1.00000000,-4.00000000):13.00000000 iterations=0 eval=3 updates=0 moved=0 personal=0 global=0 vclips=0 hits=0 draws=0 rejections=0 rng=43
• zero-budget
• status=EVAL_LIMIT best=NONE iterations=0 eval=0 updates=0 moved=0 personal=0 global=0 vclips=0 hits=0 draws=0 rejections=0 rng=43
• collapsed
• t=0 best=(4.00000000,4.00000000) f=34.00000000 current_min=34.00000000 eval=3 draws=0
• t=1 best=(4.00000000,4.00000000) f=34.00000000 current_min=34.00000000 eval=6 draws=12
• t=2 best=(4.00000000,4.00000000) f=34.00000000 current_min=34.00000000 eval=9 draws=24
• status=ITER_LIMIT best=(4.00000000,4.00000000):34.00000000 iterations=2 eval=9 updates=6 moved=0 personal=0 global=0 vclips=0 hits=0 draws=24 rejections=0 rng=1226194021

Complexity, Correctness, and Invariants

Let $k$ be the number of completed iterations. For an initialized run,

$$k=\min\left(K,\left\lfloor\frac{B-N}{N}\right\rfloor\right),\qquad E=N(k+1),\qquad U=kN,$$

where $E$ is objective evaluations and $U$ is particle updates. If $B\lt N$, no objective is evaluated. This implementation has no early convergence exit, so best, average, and worst objective landscapes use the same evaluation count at fixed $N,K,B$. Quality at that budget can differ greatly; no generic bound on iterations to a global optimum follows.

Every particle-coordinate consumes two mapped random draws. If $J$ raw draws are rejected, total raw draws are $Q=2dU+J$. For this particular generator and bound $1024$, each rejected state is followed by an accepted one: it has form $2147483647-a$ with $1\le a\le1023$, whose successor is $2147483647-48271a$, outside the rejection interval. Thus at most two raw draws are needed per mapped value. Together with $\Theta(d)$ sphere evaluation and state copying, the bounded implementation has best-, average-, and worst-case runtime $\Theta((k+1)Nd)$, including input validation.

With an external evaluator costing $C_f$ and an ordinary expected-constant-cost random source, search work is $O(NC_f+kN(d+C_f))$, plus input preparation. Population size, dimensionality, and evaluation cost all matter; increasing particles can improve coverage while reducing the number of iterations affordable under a fixed budget. A convergence-to-quality analysis requires assumptions about the objective and the exact PSO variant.

Live positions, velocities, personal memories, and the shared archive need $O(Nd)$ space. These programs retain full history, including raw/capped velocities, proposed/repaired positions, random draws, and copied memories, so actual auxiliary storage is $\Theta((k+1)Nd)$ words. Dropping traces reduces it to $O(Nd)$ plus evaluator state. Detailed output also has a cost proportional to the values printed.

The useful correctness proof is an induction over evaluated visits. Each initial memory is paired with the score computed at that position. Motion caps velocity, then repairs position, so evaluated current positions stay feasible and stored velocity magnitudes stay within their caps. Strict personal updates preserve the best visited score of each particle:

$$f(\mathbf p_i^t)=\min_{0\le s\le t}f(\mathbf x_i^s).$$

Strict swarm updates preserve the minimum across all those visits:

$$f(\mathbf g^t)=\min_i f(\mathbf p_i^t)=\min_{i,\,0\le s\le t}f(\mathbf x_i^s).$$

These equalities refer to the deterministic computed objective, with fixed arithmetic and coherent position/score pairs. Therefore personal-best and swarm-best scores never increase, and the archive's score is no greater than any current score. Current scores, the current minimum, and current mean need not improve. The archive point may be absent from the current swarm, as the worked run demonstrates.

Copying memories prevents later motion from changing their stored positions. Freezing the social guide makes each complete iteration use the specified synchronous rule. Reserving all $N$ evaluations before an iteration keeps returned state at a completed boundary. None of these invariants establishes global optimality: identical suboptimal memories and positions with zero velocities form a fixed point even when $w,c_1,c_2$ are positive.

Variants, Pitfalls, and Optimizations

Avoid memory aliasing and stale scores. Assigning a mutable current-position array directly to a personal best or swarm archive can change the remembered position during later motion without changing its cached objective. Copy the vector when recording a best, and copy nested rows when retaining a snapshot. The displayed counterexample uses the actual move from $(0,-1)$ to $(-1.5,-2.5)$: the true scores are $1$ and $8.5$, so an aliased “best” position with cached score $1$ is inconsistent.

Keep synchronous and asynchronous variants distinct. A synchronous sweep uses one old social guide for all particles. Updating and reusing that guide immediately gives an order-dependent asynchronous variant. Both can be studied, but mixing them silently changes the algorithm. For reproducible parallel evaluation, assign random draws and evaluations to fixed particle/coordinate indices, collect the full iteration, and commit memories in a defined tie-breaking order.

Specify boundary and velocity behavior together. This implementation caps velocity, proposes a position, clips strict box crossings, and zeros only the crossed velocity coordinates. Reflection, wrapping, random relocation, and alternative damping rules produce different motion. Resetting the entire velocity when only one coordinate crosses also changes the policy. Box repair alone does not satisfy nonlinear engineering constraints. A zero cap freezes that coordinate; a fixed bound may cause its initial nonzero velocity to be zeroed on the first attempted crossing.

Do not interpret attraction coefficients as probabilities. The random multipliers are in $[0,1)$, but $c_1,c_2$ can exceed one. Large pulls and momentum may cause oscillation or repeated clipping; small caps and coefficients can make progress slow. The accepted validation ranges keep arithmetic manageable here; they are not a theorem guaranteeing stable convergence for every objective or trajectory. PSO stability analyses depend on assumptions, topology, randomness, and the exact update rule. A bounded trajectory or convergent velocity is still not proof of finding a global minimizer.

Preserve diversity without claiming random attraction is exploration everywhere. If every particle, personal best, and swarm best equals the same suboptimal point and all velocities are zero, both attraction differences are zero for every random multiplier. The swarm stays there forever. Restarts, random immigrants, or hybrid local search can introduce new candidates, but they change the method and its evaluation budget. A tiny swarm spread or stalled archive is a practical stopping signal, not a certificate.

Choose the neighborhood and parameter policy deliberately. Global-best PSO shares one record with every particle. Local-best or ring variants use neighborhood records, altering information flow and exploration. Inertia schedules, coefficient adaptation, constriction formulations, and hybrid methods require their own equations and assumptions; do not add a constriction multiplier while treating it as an unrelated cosmetic parameter. The PySwarms global-best implementation is a primary reference for position/velocity updates and memory handling; its random precision, boundaries, stopping rules, and initialization policies differ from this teaching program.

Normalize and evaluate honestly. Variables with very different ranges can benefit from scaled coordinates and cap choices tied to their meaningful units. Binary, integer, and categorical decisions need appropriate PSO variants or a carefully specified decoder. Compare repeated seeds under equal objective budgets against random search and suitable local methods. A smooth convex sphere is an excellent correctness check but can be solved directly; it does not establish that a swarm is the most efficient solver.

Treat noisy objectives as a different memory problem. A lucky noisy observation can become a persistent attractor if cached indefinitely. Use a stated resampling or statistical policy and validate final candidates on independent data. Returning the archive rather than only the final current minimum preserves historical quality only when the stored point and its score remain meaningful.

Practice Problems and Real-World Use

Begin with objectives whose reference solutions are known. For the box-constrained shifted sphere, clipping the center into each coordinate interval gives the exact real-valued minimizer, since the squared terms minimize independently. Compare returned positions and computed scores with that reference, taking floating-point rounding into account. Test a nonzero initial velocity as well as a stationary swarm so the role of momentum is visible.

For an applied study, document the objective, units, bounds, constraints, initialization, coefficients, random precision and seed policy, topology, synchronous/asynchronous timing, velocity and position handling, memory ties, budget, and final validation. Report variation across seeds and objective evaluations as well as the best run. Sensor calibration and simulation-based design can supply meaningful continuous objectives; feasibility and evaluator reliability remain part of the problem.