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 Ni^{k}} \tau{il}^{\alpha}\,\eta{il}^{\beta}}, \qquad j \in Ni^{k}, $$

where tau is the learned pheromone trail and eta the static heuristic desirability (for TSP, etaij = 1/dij). 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 [taumin, taumax], restarts | robust default, best with local search |

ACS mechanics. With probability q0 the ant exploits: j = argmax over N of tauil * etail^beta; otherwise it samples with the random proportional rule. Each ant also applies a local update while constructing, tauij 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)."""

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

def traillimits(bestcost: float, rho: float, pbest: float, n: int) -> tuple[float, float]: """Return (taumin, taumax): taumax = 1/(rhoCbest), taumin from pbest.""" taumax = 1.0 / (rho bestcost) p = pbest ** (1.0 / n) taumin = taumax (1.0 - p) / ((n / 2.0 - 1.0) p) return min(taumin, taumax), tau_max

def runmmas( 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.defaultrng(cfg.seed) tau = np.full((n, n), 1.0) # rescaled to taumax after iteration 0 bestsol = np.empty(0, dtype=np.int64) bestcost = float("inf") history: list[float] = [] for it in range(cfg.niterations): itersol, itercost = bestsol, float("inf") for in range(cfg.nants): sol = construct(tau, rng) if improve is not None: sol = improve(sol) cost = evaluate(sol) if cost 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] * etab[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 = runmmas( 10, construct, evaluate, deposit, MMASConfig(nants=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]  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.

"""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 = 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.

"""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())
    f

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.