# Ant Colony Optimization

> When the user wants to design, implement, or tune ant colony optimization for combinatorial problems, including pheromone model design, visibility heuristics, and choosing among Ant System, Ant Colony System, and MAX-MIN Ant System. Also use when the user mentions "ant colony," "ACO," "pheromone," "MMAS," "ant colony system," "pheromone evaporation," "pheromone trail," "construction graph," or when solutions are built component-by-component guided by learned trail values. For tour neighborhoods and 2-opt machinery, see local-search-and-neighborhoods; for TSP formulations and benchmarks, see traveling-salesman-problem.

- Skill: `hajibabaie/ant-colony-optimization` (Agent Skill)
- Install (CLI): `npx skillmds@latest add hajibabaie/ant-colony-optimization`
- Raw SKILL.md: https://api.skillmd.com/api/skills/hajibabaie/ant-colony-optimization/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Web & Frontend
- Author: hajibabaie (https://skillmd.com/u/hajibabaie)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/hajibabaie/ant-colony-optimization

---


# Ant Colony Optimization

You are an expert in ant colony optimization (ACO) for combinatorial optimization. This skill covers
pheromone models, visibility (heuristic) information, the main variants — Ant System, Ant Colony
System (ACS), and MAX-MIN Ant System (MMAS) — pheromone trail limits, hybridization with local
search, and the design of construction graphs for problems beyond routing. Use the framework below
to select a variant, build a correct and fast implementation, and diagnose convergence behavior.

## Initial Assessment

Establish the following before writing any code or recommending a configuration:

- **Problem class and construction graph.** Identify the sequential decisions an ant makes: edge
  pheromones for routing, (item, position) pairs for assignment, per-item trails for subset
  problems. If no natural sequential construction exists, ACO is a poor fit.
- **Instance size.** The pheromone matrix is O(n^2) memory and a naive construction step is O(n).
  For n above a few thousand, candidate lists and vectorized construction are mandatory.
- **Heuristic information.** Check whether a greedy desirability measure (eta) exists, e.g. inverse
  distance for TSP. Without one (QAP, much of scheduling), plan for beta = 0 and lean on local search.
- **Local search availability.** ACO without local search is rarely competitive on classic
  benchmarks (Dorigo & Stützle 2004, Ant Colony Optimization). Confirm a delta-evaluable
  improvement procedure exists before promising results.
- **Evaluation cost.** Count objective evaluations per iteration: ants × (construction + local
  search). If the objective is expensive, shrink the colony and deepen local search.
- **Time budget and termination.** Fix a wall-clock or iteration budget up front; ACO has no
  natural stopping point. Plan stagnation detection and restarts for long runs.
- **Quality target.** Determine whether the goal is "good feasible quickly" (favor ACS, high q0) or
  "near-optimal given hours" (favor MMAS with restarts and strong local search).
- **Constraint structure.** Decide whether construction can always stay feasible (visit-once
  constraints are free in permutation construction) or whether a repair/penalty layer is needed.
- **Baseline.** Multi-start local search at the same evaluation budget is the honest baseline;
  ACO must beat it to justify its complexity.
- **Reproducibility.** Require seeded runs (`np.random.default_rng(seed)`) and multiple seeds per instance.

## Algorithm Anatomy

### Construction graph and transition rule

ACO builds solutions component by component on a construction graph. Ant k at decision point i
chooses the next component j from the feasible neighborhood N_i^k with the **random proportional
rule** of Ant System (Dorigo, Maniezzo & Colorni 1996, Ant System):

$$
p_{ij}^{k} \;=\; \frac{\tau_{ij}^{\alpha}\,\eta_{ij}^{\beta}}
{\sum_{l \in N_i^{k}} \tau_{il}^{\alpha}\,\eta_{il}^{\beta}},
\qquad j \in N_i^{k},
$$

where tau is the learned pheromone trail and eta the static heuristic desirability (for TSP,
eta_ij = 1/d_ij). Alpha and beta weight learning against greediness.

### Pheromone update

After all m ants finish, Ant System evaporates and deposits:

$$
\tau_{ij} \;\leftarrow\; (1-\rho)\,\tau_{ij} \;+\; \sum_{k=1}^{m} \Delta\tau_{ij}^{k},
\qquad
\Delta\tau_{ij}^{k} =
\begin{cases}
1/C^{k} & \text{if ant } k \text{ used component } (i,j),\\
0 & \text{otherwise,}
\end{cases}
$$

with C^k the cost of ant k's solution. Evaporation rate rho in (0, 1] controls how fast the colony
forgets. The two strong modern variants change *who deposits* and *how trails are bounded*.

### Variant comparison

| Variant | Who deposits | Selection rule | Extra mechanics | Use when |
|---|---|---|---|---|
| Ant System (1996) | all ants | random proportional | none | teaching, baselines only |
| Elitist / Rank-based AS | best + top-ranked ants | random proportional | weighted deposits | mild upgrade of AS; mostly historical |
| ACS (Dorigo & Gambardella 1997) | best-so-far only | pseudorandom proportional (q0) | local evaporation during construction | small budgets, fast convergence, online/dynamic problems |
| MMAS (Stützle & Hoos 2000) | iteration-best / best-so-far schedule | random proportional | trail limits [tau_min, tau_max], restarts | robust default, best with local search |

**ACS mechanics.** With probability q0 the ant exploits: j = argmax over N of tau_il * eta_il^beta;
otherwise it samples with the random proportional rule. Each ant also applies a **local update**
while constructing, tau_ij <- (1 - xi) tau_ij + xi tau_0 with tau_0 = 1 / (n * C_nn), which makes
edges less attractive immediately after use and decorrelates ants within one iteration. The global
update evaporates and deposits only on the best-so-far solution's components.

**MMAS mechanics.** Only the iteration-best or best-so-far ant deposits. All trails are clamped to
[tau_min, tau_max] and initialized at tau_max, which forces early exploration and provably prevents
the probability of any component from reaching 0 or 1:

$$
\tau_{\max} = \frac{1}{\rho\,C^{bs}},
\qquad
\tau_{\min} = \tau_{\max}\cdot
\frac{1 - p_{\mathrm{best}}^{1/n}}{(\,n/2 - 1\,)\;p_{\mathrm{best}}^{1/n}},
$$

where p_best (e.g. 0.05) is the target probability of reconstructing the current best solution at
full convergence, and n/2 estimates the average number of choices per step. When trails stagnate,
MMAS resets all trails to tau_max (a restart).

### Pheromone models beyond routing

The pheromone must index the decision you want the colony to learn. Standard designs:

| Problem family | Pheromone index tau | Visibility eta | Deposit on |
|---|---|---|---|
| Routing (TSP, VRP) | edge (i, j) | 1/d_ij | edges of the tour/routes |
| Assignment (QAP, timetabling) | pair (item i, location/value k) | usually none (beta = 0) | chosen (item, location) pairs |
| Subset selection (knapsack, set covering) | item i | profit/weight or cover ratio | chosen items |
| Sequencing/scheduling (flow shop, RCPSP) | (job j, position p) | priority-rule value | chosen positions, read with the summation rule |

Design rules: (1) solution quality must map monotonically to deposit (use 1/C for minimization);
(2) the trail consulted at a decision must be conditionally meaningful — in position-based
scheduling models the raw trail tau[j, p] misleads, so the summation rule (Merkle, Middendorf &
Schmeck 2002) reads sum over q <= p of tau[j, q] instead; (3) keep the pheromone O(n^2) or smaller.

### Decision guidance

- Default to **MMAS + local search on every ant**: it is the most robust published configuration
  for TSP- and QAP-like problems.
- Choose **ACS** when the evaluation budget is small, you need good solutions within seconds, or
  the problem changes online (local evaporation adapts faster).
- Per-iteration cost is O(m·n^2) for full construction scans, O(m·n·cl) with candidate lists of
  size cl, plus local search. Memory is O(n^2) for tau (and eta if dense).
- Prefer a different metaheuristic when there is no sequential construction (pure binary flips),
  when the problem is continuous, or when evaluations are too expensive for population sampling.

## Generic MMAS Framework

The only problem-specific elements of ACO are the construction procedure, the evaluation, and the
mapping from a solution to the pheromone entries it used. The framework below owns everything
else: pheromone memory, trail limits, evaporation, deposit scheduling, best-so-far bookkeeping.
For the `improve` step see **local-search-and-neighborhoods**; for cheap `evaluate` see
**fitness-evaluation-and-caching**.

```text
MMAS skeleton
-------------
initialize tau[i][j] := tau_max estimated from a quick starting solution
repeat for n_iterations:
    for each of m ants:
        s := construct(tau, rng)            # problem-specific, random proportional rule
        s := improve(s)                     # optional local search
    update iteration-best and best-so-far
    recompute tau_max := 1 / (rho * C_best),  tau_min from p_best
    tau := (1 - rho) * tau                  # evaporation
    pick deposit solution: iteration-best, every gb_every-th iteration best-so-far
    tau[components of deposit solution] += 1 / C_deposit
    clip tau into [tau_min, tau_max]; (optionally) restart tau := tau_max on stagnation
return best-so-far
```

```python
"""Reusable MAX-MIN Ant System core: pheromone memory, trail limits, updates.

Problem-specific parts are injected as callables:
  construct(tau, rng) -> solution        build one solution from trail values
  evaluate(solution)  -> float           minimization objective
  deposit(solution)   -> (rows, cols)    pheromone entries the solution used
  improve(solution)   -> solution        optional local search
"""
from __future__ import annotations

from dataclasses import dataclass
from typing import Callable

import numpy as np


@dataclass
class MMASConfig:
    """MMAS parameters (Stützle & Hoos 2000, MAX-MIN Ant System)."""

    n_ants: int = 25
    rho: float = 0.2          # evaporation rate in (0, 1]
    p_best: float = 0.05      # target probability of rebuilding the best solution
    n_iterations: int = 500
    gb_every: int = 25        # deposit best-so-far every gb_every iterations
    seed: int = 0


def trail_limits(best_cost: float, rho: float, p_best: float, n: int) -> tuple[float, float]:
    """Return (tau_min, tau_max): tau_max = 1/(rho*C_best), tau_min from p_best."""
    tau_max = 1.0 / (rho * best_cost)
    p = p_best ** (1.0 / n)
    tau_min = tau_max * (1.0 - p) / ((n / 2.0 - 1.0) * p)
    return min(tau_min, tau_max), tau_max


def run_mmas(
    n: int,
    construct: Callable[[np.ndarray, np.random.Generator], np.ndarray],
    evaluate: Callable[[np.ndarray], float],
    deposit: Callable[[np.ndarray], tuple[np.ndarray, np.ndarray]],
    cfg: MMASConfig,
    improve: Callable[[np.ndarray], np.ndarray] | None = None,
) -> tuple[np.ndarray, float, list[float]]:
    """Run MMAS; return (best solution, best cost, best-so-far history)."""
    rng = np.random.default_rng(cfg.seed)
    tau = np.full((n, n), 1.0)            # rescaled to tau_max after iteration 0
    best_sol = np.empty(0, dtype=np.int64)
    best_cost = float("inf")
    history: list[float] = []
    for it in range(cfg.n_iterations):
        iter_sol, iter_cost = best_sol, float("inf")
        for _ in range(cfg.n_ants):
            sol = construct(tau, rng)
            if improve is not None:
                sol = improve(sol)
            cost = evaluate(sol)
            if cost < iter_cost:
                iter_sol, iter_cost = sol, cost
        if iter_cost < best_cost:
            best_sol, best_cost = iter_sol.copy(), iter_cost
        tau_min, tau_max = trail_limits(best_cost, cfg.rho, cfg.p_best, n)
        if it == 0:
            tau.fill(tau_max)             # MMAS initialization at the upper limit
        tau *= 1.0 - cfg.rho
        use_gb = (it + 1) % cfg.gb_every == 0
        dep_sol, dep_cost = (best_sol, best_cost) if use_gb else (iter_sol, iter_cost)
        rows, cols = deposit(dep_sol)
        tau[rows, cols] += 1.0 / dep_cost
        np.clip(tau, tau_min, tau_max, out=tau)
        history.append(best_cost)
    return best_sol, best_cost, history


if __name__ == "__main__":
    # Tiny demo: 10-city Euclidean TSP wired into the framework via closures.
    rng0 = np.random.default_rng(42)
    pts = rng0.random((10, 2))
    dist = np.sqrt(((pts[:, None, :] - pts[None, :, :]) ** 2).sum(axis=2))
    eta_b = (1.0 / np.maximum(dist, 1e-9)) ** 2.0          # beta = 2 folded into eta

    def construct(tau: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        n = tau.shape[0]
        tour = np.empty(n, dtype=np.int64)
        visited = np.zeros(n, dtype=bool)
        cur = int(rng.integers(n))
        tour[0] = cur
        visited[cur] = True
        for step in range(1, n):
            cand = np.flatnonzero(~visited)
            w = tau[cur, cand] * eta_b[cur, cand]          # alpha = 1
            cur = int(cand[rng.choice(cand.size, p=w / w.sum())])
            tour[step] = cur
            visited[cur] = True
        return tour

    def evaluate(tour: np.ndarray) -> float:
        return float(dist[tour, np.roll(tour, -1)].sum())

    def deposit(tour: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        nxt = np.roll(tour, -1)                            # symmetric deposit
        return np.concatenate([tour, nxt]), np.concatenate([nxt, tour])

    sol, cost, hist = run_mmas(
        10, construct, evaluate, deposit, MMASConfig(n_ants=10, n_iterations=120, seed=7)
    )
    print(f"10-city tour length: {cost:.4f}")
    # Expected: a valid 10-city tour; length typically 2.5-3.2 and equal to the
    # optimum for most seeds (10 uniform points; exhaustive check is feasible).
```

### Parameter guidance

| Parameter | Typical range | What it trades off |
|---|---|---|
| m (ants per iteration) | 10–50; with local search 10–25 | sampling breadth per iteration vs. iterations within the budget |
| alpha (trail weight) | 1.0, rarely tuned | >1 amplifies trail differences — faster convergence, higher stagnation risk |
| beta (visibility weight) | 2–5 (0 when no eta exists) | high beta = greedy and strong early; low beta lets learned trails dominate |
| rho (evaporation) | MMAS: 0.2 with local search, 0.02 without; ACS: 0.1 | fast forgetting adapts quickly but is unstable; slow forgetting learns slowly |
| q0 (ACS exploitation) | 0.7–0.95 | exploitation of the best-known components vs. exploration during construction |
| xi (ACS local update) | 0.1 | stronger local evaporation decorrelates ants inside one iteration |
| p_best (MMAS) | 0.005–0.05 | smaller value raises tau_min — more exploration, slower convergence |
| cl (candidate list size) | 10–30 | smaller is faster and usually better; too small can exclude needed edges |
| gb_every (deposit schedule) | start 25, decrease over the run | iteration-best favors exploration; best-so-far accelerates convergence |

## Worked Example 1: TSP with MMAS + 2-opt

The reference configuration from Stützle & Hoos (2000): MMAS with a nearest-neighbor candidate
list, 2-opt applied to every ant's tour, deposit alternating between iteration-best and
best-so-far. The 2-opt here is a compact best-improvement version with a vectorized move scan; for
don't-look bits, neighbor-list pruning, and Or-opt see **local-search-and-neighborhoods**.

```python
"""TSP utilities for the MMAS worked example: instance, 2-opt, candidate lists."""
from __future__ import annotations

import numpy as np


def euclidean_instance(n: int, seed: int) -> np.ndarray:
    """n random points in the unit square -> (n, n) Euclidean distance matrix."""
    rng = np.random.default_rng(seed)
    pts = rng.random((n, 2))
    diff = pts[:, None, :] - pts[None, :, :]
    return np.sqrt((diff * diff).sum(axis=2))


def tour_length(tour: np.ndarray, dist: np.ndarray) -> float:
    """Cyclic tour length."""
    return float(dist[tour, np.roll(tour, -1)].sum())


def two_opt(tour: np.ndarray, dist: np.ndarray) -> np.ndarray:
    """Best-improvement 2-opt with a vectorized move scan per anchor edge."""
    tour = tour.copy()
    n = tour.size
    improved = True
    while improved:
        improved = False
        succ = np.roll(tour, -1)
        d_cur = dist[tour, succ]
        for i in range(n - 2):
            a, b = tour[i], tour[i + 1]
            j_hi = n - 1 if i == 0 else n      # skip the move that flips the whole tour
            js = np.arange(i + 2, j_hi)
            if js.size == 0:
                continue
            c, d = tour[js], succ[js]
            delta = dist[a, c] + dist[b, d] - d_cur[i] - d_cur[js]
            k = int(np.argmin(delta))
            if delta[k] < -1e-12:
                j = int(js[k])
                tour[i + 1 : j + 1] = tour[i + 1 : j + 1][::-1]
                improved = True
                succ = np.roll(tour, -1)
                d_cur = dist[tour, succ]
    return tour


def nn_candidate_lists(dist: np.ndarray, cl: int) -> np.ndarray:
    """cl nearest neighbors per city, excluding the city itself; shape (n, cl)."""
    order = np.argsort(dist, axis=1)
    return order[:, 1 : cl + 1]


if __name__ == "__main__":
    dist = euclidean_instance(30, seed=0)
    rng = np.random.default_rng(1)
    rand_tour = rng.permutation(30).astype(np.int64)
    opt_tour = two_opt(rand_tour, dist)
    print(f"random: {tour_length(rand_tour, dist):.3f} -> 2-opt: {tour_length(opt_tour, dist):.3f}")
    # Expected: 2-opt shortens a random tour substantially, e.g. ~14-17 down to ~4.4-5.0.
```

The MMAS driver below continues the same module (it calls the helpers defined above; imports are
repeated so the block parses on its own). The weight matrix `tau**alpha * eta**beta` is recomputed
once per iteration and shared by all ants — trails do not change during MMAS construction.

```python
"""MMAS for the symmetric TSP with 2-opt on every ant (same module as above)."""
from __future__ import annotations

import numpy as np


def construct_tour(w: np.ndarray, cand: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """One ant: random proportional rule on the candidate list, full-scan fallback."""
    n = w.shape[0]
    tour = np.empty(n, dtype=np.int64)
    visited = np.zeros(n, dtype=bool)
    cur = int(rng.integers(n))
    tour[0] = cur
    visited[cur] = True
    for step in range(1, n):
        nbrs = cand[cur]
        nbrs = nbrs[~visited[nbrs]]
        if nbrs.size == 0:                       # candidate list exhausted
            nbrs = np.flatnonzero(~visited)
        weights = w[cur, nbrs]
        total = weights.sum()
        if total > 0.0:
            cur = int(nbrs[rng.choice(nbrs.size, p=weights / total)])
        else:
            cur = int(nbrs[rng.integers(nbrs.size)])
        tour[step] = cur
        visited[cur] = True
    return tour


def mmas_tsp(
    dist: np.ndarray,
    n_ants: int = 20,
    alpha: float = 1.0,
    beta: float = 3.0,
    rho: float = 0.2,
    p_best: float = 0.05,
    cl: int = 15,
    n_iterations: int = 300,
    gb_every: int = 25,
    seed: int = 0,
) -> tuple[np.ndarray, float, list[float]]:
    """MAX-MIN Ant System for symmetric TSP, 2-opt applied to every ant's tour."""
    n = dist.shape[0]
    rng = np.random.default_rng(seed)
    eta_beta = (1.0 / np.maximum(dist, 1e-12)) ** beta
    np.fill_diagonal(eta_beta, 0.0)
    cand = nn_candidate_lists(dist, min(cl, n - 1))
    start = two_opt(rng.permutation(n).astype(np.int64), dist)
    best_tour, best_cost = start, tour_length(start, dist)
    tau = np.full((n, n), 1.0 / (rho * best_cost))         # tau_max initialization
    history: list[float] = []
    for it in range(n_iterations):
        w = (tau ** alpha) * eta_beta                      # shared by all ants
        iter_tour, iter_cost = best_tour, float("inf")
        for _ in range(n_ants):
            tour = two_opt(construct_tour(w, cand, rng), dist)
            cost = tour_length(tour, dist)
            if cost < iter_cost:
                iter_tour, iter_cost = tour, cost
        if iter_cost < best_cost:
            best_tour, best_cost = iter_tour.copy(), iter_cost
        tau_max = 1.0 / (rho * best_cost)
        p = p_best ** (1.0 / n)
        tau_min = min(tau_max, tau_max * (1.0 - p) / ((n / 2.0 - 1.0) * p))
        tau *= 1.0 - rho
        dep = best_tour if (it + 1) % gb_every == 0 else iter_tour
        dep_cost = tour_length(dep, dist)
        nxt = np.roll(dep, -1)
        tau[dep, nxt] += 1.0 / dep_cost                    # symmetric deposit
        tau[nxt, dep] += 1.0 / dep_cost
        np.clip(tau, tau_min, tau_max, out=tau)
        history.append(best_cost)
    return best_tour, best_cost, history


if __name__ == "__main__":
    dist = euclidean_instance(40, seed=42)
    tour, cost, hist = mmas_tsp(dist, n_iterations=120, seed=1)
    assert sorted(tour.tolist()) == list(range(40))
    assert all(a >= b for a, b in zip(hist, hist[1:]))     # best-so-far never worsens
    print(f"best tour length: {cost:.4f}")
    # Expected: a valid 40-city tour with length ~5.0-5.5; the best-so-far history
    # is non-increasing and the final value matches or beats multi-start 2-opt.
```

Implementation notes: candidate lists change construction from O(n) to O(cl) per step and usually
*improve* quality, because long edges are excluded from sampling. Initialization at tau_max keeps
early iterations close to a randomized greedy heuristic; learning takes over as trails
differentiate. For symmetric TSP deposit on both (i, j) and (j, i); for asymmetric TSP deposit
only on the directed arcs actually traversed.

## Variant Mechanics: Ant Colony System

ACS converges faster than MMAS at the price of less exploration; it is the better choice for
tight time budgets. The block below isolates ACS without local search so the three differences are
visible: pseudorandom proportional rule, in-construction local update, best-so-far-only global update.

```python
"""Ant Colony System for symmetric TSP (Dorigo & Gambardella 1997)."""
from __future__ import annotations

import numpy as np


def nearest_neighbor_cost(dist: np.ndarray) -> float:
    """Length of the nearest-neighbor tour from city 0, used to calibrate tau_0."""
    n = dist.shape[0]
    visited = np.zeros(n, dtype=bool)
    visited[0] = True
    cur, total = 0, 0.0
    for _ in range(n - 1):
        d = np.where(visited, np.inf, dist[cur])
        nxt = int(np.argmin(d))
        total += dist[cur, nxt]
        visited[nxt] = True
        cur = nxt
    return total + dist[cur, 0]


def acs_tsp(
    dist: np.ndarray,
    n_ants: int = 10,
    beta: float = 2.0,
    rho: float = 0.1,
    xi: float = 0.1,
    q0: float = 0.9,
    n_iterations: int = 400,
    seed: int = 0,
) -> tuple[np.ndarray, float]:
    """ACS without local search, to isolate the variant's mechanics."""
    n = dist.shape[0]
    rng = np.random.default_rng(seed)
    eta_beta = (1.0 / np.maximum(dist, 1e-12)) ** beta
    np.fill_diagonal(eta_beta, 0.0)
    tau0 = 1.0 / (n * nearest_neighbor_cost(dist))
    tau = np.full((n, n), tau0)
    best_tour = np.arange(n, dtype=np.int64)
    best_cost = float(dist[best_tour, np.roll(best_tour, -1)].sum())
    for _ in range(n_iterations):
        for _ in range(n_ants):
            tour = np.empty(n, dtype=np.int64)
            visited = np.zeros(n, dtype=bool)
            cur = int(rng.integers(n))
            tour[0] = cur
            visited[cur] = True
            for step in range(1, n):
                cand = np.flatnonzero(~visited)
                scores = tau[cur, cand] * eta_beta[cur, cand]
                if rng.random() < q0:                      # exploit
                    nxt = int(cand[np.argmax(scores)])
                else:                                      # biased exploration
                    total = scores.sum()
                    if total > 0.0:
                        nxt = int(cand[rng.choice(cand.size, p=scores / total)])
                    else:
                        nxt = int(rng.choice(cand))
                tau[cur, nxt] = (1.0 - xi) * tau[cur, nxt] + xi * tau0   # local update
                tau[nxt, cur] = tau[cur, nxt]
                tour[step] = nxt
                visited[nxt] = True
                cur = nxt
            cost = float(dist[tour, np.roll(tour, -1)].sum())
            if cost < best_cost:
                best_tour, best_cost = tour.copy(), cost
        nxt_b = np.roll(best_tour, -1)                     # global update: best edges only
        tau[best_tour, nxt_b] = (1.0 - rho) * tau[best_tour, nxt_b] + rho / best_cost
        tau[nxt_b, best_tour] = tau[best_tour, nxt_b]
    return best_tour, best_cost


if __name__ == "__main__":
    rng = np.random.default_rng(5)
    pts = rng.random((30, 2))
    dist = np.sqrt(((pts[:, None, :] - pts[None, :, :]) ** 2).sum(axis=2))
    tour, cost = acs_tsp(dist, n_iterations=250, seed=2)
    assert sorted(tour.tolist()) == list(range(30))
    print(f"ACS tour length: {cost:.4f}")
    # Expected: a valid 30-city tour, length ~4.4-5.0 without local search; adding
    # 2-opt per ant (as in the MMAS example) closes most of the remaining gap.
```

Note the asymmetry of the global update: in ACS, evaporation touches only the best-so-far edges.
This is what makes ACS aggressive — the matrix drifts toward tau0 while elite edges are reinforced.

## Worked Example 2: QAP with Pair-Assignment Pheromones

The quadratic assignment problem places n facilities at n locations to minimize the flow-weighted
distance sum: min over permutations p of sum_{i,j} flow[i,j] * dist[p[i], p[j]]. The pheromone here
is **not** on edges of a path: tau[i, k] is the learned desirability of placing facility i at
location k — the pair-assignment model. QAP has no reliable a-priori visibility, so beta = 0,
following MMAS-QAP (Stützle & Hoos 2000): pheromone-only construction plus 2-exchange local search
with O(n) delta evaluation (see **fitness-evaluation-and-caching** for the general pattern).

```python
"""QAP utilities: instance generator, vectorized cost, exact O(n) swap delta."""
from __future__ import annotations

import numpy as np


def random_qap_instance(n: int, seed: int) -> tuple[np.ndarray, np.ndarray]:
    """Symmetric random QAP, integer flows/distances in [0, 100), zero diagonal."""
    rng = np.random.default_rng(seed)
    flow = np.triu(rng.integers(0, 100, (n, n)), 1)
    dist = np.triu(rng.integers(0, 100, (n, n)), 1)
    return (flow + flow.T).astype(np.float64), (dist + dist.T).astype(np.float64)


def qap_cost(p: np.ndarray, flow: np.ndarray, dist: np.ndarray) -> float:
    """sum_{i,j} flow[i,j] * dist[p[i], p[j]], fully vectorized."""
    return float((flow * dist[np.ix_(p, p)]).sum())


def swap_delta(p: np.ndarray, flow: np.ndarray, dist: np.ndarray, r: int, s: int) -> float:
    """Exact O(n) cost change of swapping the locations of facilities r and s.

    Valid for asymmetric flow/distance matrices as well (Taillard 1991 formula).
    """
    pr, ps = p[r], p[s]
    mask = np.ones(p.size, dtype=bool)
    mask[[r, s]] = False
    pk = p[mask]
    delta = (
        (flow[r, r] - flow[s, s]) * (dist[ps, ps] - dist[pr, pr])
        + (flow[r, s] - flow[s, r]) * (dist[ps, pr] - dist[pr, ps])
        + (flow[r, mask] - flow[s, mask]) @ (dist[ps, pk] - dist[pr, pk])
        + (flow[mask, r] - flow[mask, s]) @ (dist[pk, ps] - dist[pk, pr])
    )
    return float(delta)


def two_exchange(p: np.ndarray, flow: np.ndarray, dist: np.ndarray) -> np.ndarray:
    """First-improvement 2-exchange local search; loops until locally optimal."""
    p = p.copy()
    n = p.size
    improved = True
    while improved:
        improved = False
        for r in range(n - 1):
            for s in range(r + 1, n):
                if swap_delta(p, flow, dist, r, s) < -1e-9:
                    p[r], p[s] = p[s], p[r]
                    improved = True
    return p


if __name__ == "__main__":
    flow, dist = random_qap_instance(8, seed=1)
    rng = np.random.default_rng(2)
    p = rng.permutation(8).astype(np.int64)
    base = qap_cost(p, flow, dist)
    q = p.copy()
    q[2], q[5] = q[5], q[2]
    assert abs((qap_cost(q, flow, dist) - base) - swap_delta(p, flow, dist, 2, 5)) < 1e-9
    print("delta evaluation matches full recomputation")
    # Expected: the assertion passes -- the O(n) delta equals the O(n^2) recomputation.
```

The ACO driver continues the same module. Construction assigns facilities in a *random order* each
time — a fixed order would bias early facilities toward locations chosen under an emptier board —
and samples a free location for facility i in proportion to tau[i, :].

```python
"""MMAS for the QAP with pair-assignment pheromones (same module as above)."""
from __future__ import annotations

import numpy as np


def construct_assignment(tau: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Assign facilities in random order; sample a free location ~ tau[i, :]."""
    n = tau.shape[0]
    p = np.empty(n, dtype=np.int64)
    free = np.ones(n, dtype=bool)
    for i in rng.permutation(n):
        w = np.where(free, tau[i], 0.0)
        total = w.sum()
        if total > 0.0:
            loc = int(rng.choice(n, p=w / total))
        else:
            loc = int(rng.choice(np.flatnonzero(free)))
        p[i] = loc
        free[loc] = False
    return p


def mmas_qap(
    flow: np.ndarray,
    dist: np.ndarray,
    n_ants: int = 10,
    rho: float = 0.2,
    p_best: float = 0.05,
    n_iterations: int = 200,
    gb_every: int = 20,
    seed: int = 0,
) -> tuple[np.ndarray, float, list[float]]:
    """MMAS-QAP: beta = 0 (no visibility), 2-exchange local search on every ant."""
    n = flow.shape[0]
    rng = np.random.default_rng(seed)
    start = two_exchange(rng.permutation(n).astype(np.int64), flow, dist)
    best_p, best_cost = start, qap_cost(start, flow, dist)
    tau = np.full((n, n), 1.0 / (rho * best_cost))         # tau_max initialization
    fac = np.arange(n)
    history: list[float] = []
    for it in range(n_iterations):
        iter_p, iter_cost = best_p, float("inf")
        for _ in range(n_ants):
            p = two_exchange(construct_assignment(tau, rng), flow, dist)
            cost = qap_cost(p, flow, dist)
            if cost < iter_cost:
                iter_p, iter_cost = p, cost
        if iter_cost < best_cost:
            best_p, best_cost = iter_p.copy(), iter_cost
        tau_max = 1.0 / (rho * best_cost)
        q = p_best ** (1.0 / n)
        tau_min = min(tau_max, tau_max * (1.0 - q) / ((n / 2.0 - 1.0) * q))
        tau *= 1.0 - rho
        dep = best_p if (it + 1) % gb_every == 0 else iter_p
        tau[fac, dep] += 1.0 / qap_cost(dep, flow, dist)   # deposit on (i, p[i]) pairs
        np.clip(tau, tau_min, tau_max, out=tau)
        history.append(best_cost)
    return best_p, best_cost, history


if __name__ == "__main__":
    flow, dist = random_qap_instance(12, seed=7)
    p, cost, hist = mmas_qap(flow, dist, n_iterations=80, seed=3)
    assert sorted(p.tolist()) == list(range(12))
    print(f"best QAP cost: {cost:.0f}")
    # Expected: a valid permutation of 12 locations; cost matches the best value
    # found by ~100 independent 2-exchange restarts, usually within 30 iterations.
```

Differences from the TSP model worth internalizing: the pheromone row index is the decision
stage's identity (facility), not the previous choice — there is no "current node" and the
construction graph is bipartite. Deposit hits exactly n entries (one per facility), not n edges of
a cycle. And with beta = 0 the colony learns purely from solution quality; local search supplies
the intensification that eta supplies in routing problems.

## Advanced Techniques

### Stagnation detection and restarts

The standard convergence diagnostic is the average lambda-branching factor (Gambardella & Dorigo
1995, Ant-Q): per pheromone row, count entries above a threshold interpolated between the row's
min and max; a value near 1 means construction has become deterministic. Restart only when the
branching factor is low *and* the best-so-far has not improved for a fixed number of iterations.

```python
"""Stagnation diagnosis via the average lambda-branching factor."""
from __future__ import annotations

import numpy as np


def avg_branching_factor(tau: np.ndarray, lam: float = 0.05) -> float:
    """Mean number of 'attractive' choices per construction step (>=1)."""
    t_min = tau.min(axis=1, keepdims=True)
    t_max = tau.max(axis=1, keepdims=True)
    thresh = t_min + lam * (t_max - t_min)
    return float((tau >= thresh).sum(axis=1).mean())


def maybe_restart(
    tau: np.ndarray, tau_max: float, stale_iters: int,
    lam: float = 0.05, bf_limit: float = 2.0, stale_limit: int = 100,
) -> bool:
    """Reset all trails to tau_max when converged AND stuck; return True if reset."""
    if stale_iters >= stale_limit and avg_branching_factor(tau, lam) < bf_limit:
        tau.fill(tau_max)
        return True
    return False
    # Usage: call every iteration with stale_iters = iterations since last
    # best-so-far improvement; after a restart, deposit only iteration-best for
    # a while so the colony re-explores instead of re-converging instantly.
```

### Vectorized batch construction

For large colonies the Python-level loop over ants dominates. Advance all m ants one step at a
time and sample m categorical distributions in a single shot with the Gumbel-max trick. This needs
strictly positive weights for all unvisited components (true whenever tau >= tau_min > 0).

```python
"""Batch construction: all m ants advance one step per Python-level iteration."""
from __future__ import annotations

import numpy as np


def construct_tours_batch(w: np.ndarray, m: int, rng: np.random.Generator) -> np.ndarray:
    """Build m tours at once from weight matrix w = tau**alpha * eta**beta."""
    n = w.shape[0]
    tours = np.empty((m, n), dtype=np.int64)
    visited = np.zeros((m, n), dtype=bool)
    cur = rng.integers(n, size=m)
    tours[:, 0] = cur
    visited[np.arange(m), cur] = True
    for step in range(1, n):
        scores = np.where(visited, 0.0, w[cur])            # (m, n) gather + mask
        gumbel = -np.log(-np.log(rng.random((m, n)) + 1e-300) + 1e-300)
        with np.errstate(divide="ignore"):
            logits = np.log(scores)                        # log(0) -> -inf is intended
        cur = np.argmax(logits + gumbel, axis=1)           # categorical sample per row
        tours[:, step] = cur
        visited[np.arange(m), cur] = True
    return tours


if __name__ == "__main__":
    rng = np.random.default_rng(0)
    w = rng.random((50, 50)) + 0.1
    np.fill_diagonal(w, 0.0)
    tours = construct_tours_batch(w, m=20, rng=rng)
    assert all(sorted(t) == list(range(50)) for t in tours.tolist())
    print("20 valid tours built in", tours.shape[1] - 1, "batched steps")
    # Expected: all 20 tours are permutations of 0..49; runtime scales with n,
    # not with m*n, at the Python level.
```

### Local search hybridization schedules

Three budget knobs control the ACO/local-search balance: *which ants* get improved (all ants is
the published default for MMAS; improving only the iteration-best saves 90%+ of the local-search
time at a quality cost), *how deep* (2-opt to local optimality vs. a capped number of passes), and
*which neighborhood* (2-opt for tours, 2-exchange for assignments; 3-opt or Or-opt only on the
best-so-far solution). A robust scheme for expensive objectives: a cheap first-improvement pass
for every ant, full optimization only for the iteration-best. Always deposit the *improved*
solutions — depositing pre-local-search solutions discards what local search just learned.

### The hyper-cube framework

Blum & Dorigo (2004, the hyper-cube framework for ant colony optimization) normalize the deposit
so trails always stay in [0, 1]: tau <- (1 - rho) tau + rho * delta, where delta is 1 on the
components of the deposit solution and 0 elsewhere (or a weighted average over several solutions).
The pheromone scale becomes independent of objective magnitude, so one parameter setting transfers
across instances and numerical under/overflow disappears.

### Choosing the deposit schedule

MMAS quality is sensitive to the iteration-best vs. best-so-far mix. The published long-run
schedule tightens over time: deposit best-so-far every 25 iterations early, every 5 mid-run, every
iteration late; after a restart, return to iteration-best only. If you tune one ACO-specific
parameter beyond rho and beta, tune this schedule.

## Practical Challenges

**All ants build the same solution after a few hundred iterations.** This is trail stagnation, the
canonical ACO failure. Verify trail limits are actually applied (a missing `np.clip` is the most
common bug), lower p_best to raise tau_min, slow the deposit schedule (more iteration-best), and
add branching-factor-triggered restarts. If stagnation happens within tens of iterations, rho is
too high or alpha > 1.

**ACO is slower than multi-start local search at equal quality.** Expected when the pheromone adds
nothing — check that the best-so-far keeps improving after the first restart-free phase. If not,
the pheromone model may not match the decision structure (e.g., edge pheromones on a problem where
position matters), or eta dominates (beta too high) so the colony never learns. ACO earns its cost
when trails concentrate search on shared building blocks of elite solutions.

**Construction dominates runtime in Python.** Precompute the weight matrix tau^alpha * eta^beta
once per iteration (MMAS), use candidate lists, and switch to batched construction (Gumbel-max) or
numba for the inner loop. Profile first: with strong local search, 2-opt usually dwarfs
construction, and the right fix is don't-look bits, not faster sampling.

**Deposit magnitudes differ wildly across instances.** With deposit 1/C, an instance with costs in
the millions learns trails 10^6 times smaller than a toy instance, breaking tuned tau limits.
MMAS's coupled limits (both scale with 1/C_best) handle this; otherwise normalize costs or adopt
the hyper-cube framework.

**Asymmetric problems with a symmetric update.** For asymmetric TSP or directed routing, deposit
only on traversed arcs and keep tau asymmetric. Symmetrizing the deposit on a directed problem
blurs the learning signal; the colony learns "this edge is good in some direction," which is much
weaker information.

**The heuristic term is useless or misleading.** For QAP, timetabling, and many scheduling
problems no good static eta exists. Set beta = 0 rather than inventing a weak eta — a misleading
visibility term actively fights the learned trails. Compensate with local search on every ant.

**Trails saturate at the bounds and learning stops.** If most of tau sits at tau_min with a few
entries at tau_max, the effective search is a randomized greedy around one solution. Decrease rho
(slower convergence), increase p_best slightly, or shorten the gb_every schedule less aggressively.
Log the fraction of entries at each bound per iteration; it is the cheapest health metric.

## Tools & Libraries

| Library / tool | When to use | Note |
|---|---|---|
| numpy | always | trail matrices, vectorized weights, batch sampling |
| scipy.spatial.distance (cdist) | building distance matrices | faster and cleaner than manual broadcasting for large n |
| numba | construction loop too slow, n > ~2000 | JIT the per-ant loop; keep numpy for updates |
| ACOTSP (Stützle, C code) | reference results on TSPLIB | the canonical implementation behind the MMAS papers |
| scikit-opt (`sko.ACA`) | quick prototype on small TSPs | basic Ant System only; no trail limits or restarts |
| MEALPY | comparing many metaheuristics in one harness | ACO variants included; verify update rules before publishing numbers |
| OR-Tools routing | you need good VRP/TSP solutions, not ACO research | uses guided local search, not ACO — often the pragmatic choice |
| tsplib95 | loading TSPLIB instances | pairs with QAPLIB files (plain text) for QAP experiments |

## Output Format

A complete ACO deliverable contains:

1. **Configuration summary** — every parameter, explicit and seeded:

| Field | Example |
|---|---|
| variant | MMAS |
| ants / alpha / beta / rho | 20 / 1.0 / 3.0 / 0.2 |
| trail limits | p_best = 0.05, coupled to best-so-far |
| deposit schedule | iteration-best; best-so-far every 25 |
| local search | 2-opt, all ants, first-improvement; 15-NN candidate list |
| budget / seeds | 300 iterations; seeds 0-9 |

2. **Solution report** — best objective per seed, validated by an independent checker (re-compute
   the objective from the raw solution; assert the solution is a permutation / covers all
   customers), plus gap to best-known or to an exact bound where available.
3. **Convergence summary** — best-so-far curve data per seed (iteration, value), itera

…(truncated)