Tabu Search

Tabu Search

You are an expert in metaheuristic optimization, specifically in designing and implementing tabu search (Glover 1986, "Future paths for integer programming and links to artificial intelligence"; Glover & Laguna 1997, Tabu Search). This skill covers the memory structures that make tabu search work — tabu lists and tenure, move attributes, aspiration criteria, frequency-based long-term memory, candidate list strategies, and elite-based intensification — plus two reference implementations: Taillard's robust tabu search for the QAP and a critical-path tabu search for job-shop scheduling. Use the framework below to pick each memory component deliberately, implement against the verified code, and validate results against exact baselines on small instances.

Initial Assessment

Establish the following before writing any tabu search code. Each answer changes a design decision downstream.

• Problem class and representation. Permutation, binary vector, assignment, sequencing per machine? The representation fixes which neighborhoods and which move attributes are available. See solution-representation guidance in local-search-and-neighborhoods before adding memory on top.
• Neighborhood and its size. Tabu search scans many moves per iteration (often the whole neighborhood). An O(n²) swap neighborhood with O(1) delta lookup is the sweet spot; an O(n²) neighborhood with O(n) evaluation each already costs O(n³) per iteration and needs candidate lists.
• Delta evaluation cost. Tabu search is only competitive when move evaluation is incremental. If each move needs a full objective recomputation, fix that first (fitness-evaluation-and-caching) — memory structures cannot compensate for a slow scan.
• What should be forbidden? Decide the move attribute: the exact inverse move, an element-position assignment, a destroyed arc/edge. This is the single most important design choice (see the granularity table below).
• Tenure regime. Fixed, randomized per move, or reactive? Randomized tenure (Taillard 1991) is the robust default; fixed tenure needs per-instance tuning.
• Evaluation budget and time limit. Iterations × (scan cost) must fit the budget. Tabu search has no natural stopping point; pick max iterations, max iterations without improvement, or wall-clock time explicitly.
• Deterministic or stochastic runs? Tabu search with fixed tenure and deterministic tie-breaking is deterministic given the start. With randomized tenure/tie-breaking, plan a multi-seed protocol for reporting.
• Feasibility handling. Will the search stay inside the feasible region (feasibility-preserving moves), or oscillate across the boundary with penalties? Strategic oscillation needs a penalty schedule.
• Known optima or best-known values. For QAPLIB, Taillard, OR-Library instances, best-known values exist — report gaps against them. For new problems, build a small-instance exact baseline (MIP/CP) first.
• Cycling risk indicators. Plateaus and symmetric solutions raise cycling risk; plan solution hashing for cycle detection and a diversification mechanism from the start, not as an afterthought.
• Single long run vs restarts. Decide whether the budget goes into one long run with diversification phases or several restarts from elite solutions; both need the elite pool machinery below.
• Comparison baselines. At minimum: best-improvement local search with random restarts, and one alternative metaheuristic (e.g., simulated annealing). A tabu list must demonstrably beat memoryless descent on your instances.

Algorithm Anatomy

Tabu search is best-improvement local search plus memory. At iteration $t$, with current solution $x_t$, neighborhood $\mathcal{N}(x_t)$, active tabu attribute set $T_t$, and incumbent $x^*$:

$$ x_{t+1} = \operatorname*{argmin}_{y \in \mathcal{N}(x_t)} \left{ f(y) ;:; \text{attr}(x_t \to y) \notin T_t ;\text{ or }; f(y) < f(x^*) \right} $$

The move is executed even if $f(x_{t+1}) > f(x_t)$ — tabu search climbs out of local optima deterministically, unlike simulated annealing's probabilistic acceptance. The tabu list then forbids reversing recent moves for tenure $\tau$ iterations, so the search cannot fall straight back into the optimum it just left. The condition $f(y) < f(x^*)$ is the standard improved-best aspiration: a tabu status is overridden whenever the move reaches a new global best, because the tabu list exists to prevent cycling, and a strictly better solution cannot be a repeat.

The three memories

A minimal tabu search uses only recency memory. A production tabu search uses all three: recency every iteration, frequency when the search stagnates, quality when the budget allows restarts. Glover & Laguna (1997) treat this layering as the defining feature of the method.

Move attributes: what exactly becomes tabu

The tabu list stores attributes of executed moves, not full solutions. The attribute granularity sets how aggressively the list prunes:

Rule of thumb: pick the attribute that matches what the objective depends on (edges for tours, assignments for QAP, arcs for schedules). Too-fine attributes let the search cycle through near-duplicates; too-coarse attributes forbid large unexplored regions and force constant aspiration overrides.

Tenure

Tenure $\tau$ is how long an attribute stays tabu, in iterations. The classic guidance:

• Fixed: $\tau = 7$ was Glover's early constant; $\tau \approx \sqrt{|\mathcal{N}|}$ or $\tau \in [0.5\sqrt{n}, 2\sqrt{n}]$ are common starting points. Fixed tenure is fragile — one value rarely fits all instances.
• Randomized per move: draw $\tau \sim U[\tau_{\min}, \tau_{\max}]$ each time an attribute is recorded. Taillard (1991, "Robust taboo search for the quadratic assignment problem") uses $U[0.9n,,1.1n]$ for QAP and shows the randomization itself prevents cycles that any fixed value permits. This is the recommended default.
• Reactive: adapt $\tau$ from observed solution repetitions (Battiti & Tecchiolli 1994, "The reactive tabu search") — see Advanced Techniques.

Symptoms: tenure too short → the search returns to the same local optimum every few iterations (detect by hashing solutions); tenure too long → objective trace shows long flat ramps because most of the neighborhood is forbidden, and the all-tabu fallback fires often.

Aspiration criteria

1. Improved-best (always implement): override tabu if the move yields $f < f(x^*)$.
2. Default / least-tabu (always implement): if every move is tabu and none aspirates, execute the move whose tabu status expires earliest (or the best tabu move). Never skip an iteration — a frozen search wastes the remaining budget.
3. Forced long-term aspiration (Ro-TS): execute a move regardless of its value if it creates an assignment not seen for $A$ iterations (e.g., $A = 5n^2$). This is diversification built into the aspiration test.
4. Regional aspiration: override tabu if the move improves on the best solution within the current search region (requires defining regions; rarely worth it before the first three).

When to choose tabu search

Per-iteration complexity = (number of candidate moves scanned) × (delta lookup cost) + (memory update cost). The QAP implementation below achieves O(n²) per iteration: O(n²) scan with O(1) lookups plus O(n²) delta-matrix maintenance.

Generic Tabu Search Framework

The skeleton every variant shares:

TABU SEARCH (minimization)
Input: initial solution x0, neighborhood N(.), tenure rule tau(.),
       aspiration test A(.), iteration budget T
x <- x0;  x_best <- x0
tabu <- empty map {attribute -> iteration until which it is active}
for t = 1 .. T:
    C <- candidate moves from N(x)              # full scan or candidate list
    admissible <- { m in C : tabu[attr_check(m)] < t  or  A(m) holds }
    if admissible is empty:
        m* <- move in C whose tabu status expires earliest   # default aspiration
    else:
        m* <- argmin_{m in admissible} f(x (+) m)   # best admissible, even uphill
    x <- x (+) m*
    tabu[attr_record(m*)] <- t + tau(t)         # forbid reversal for tau iterations
    if f(x) < f(x_best): x_best <- x
    # optional: update frequency memory; trigger diversification/intensification
return x_best

The generic engine below separates the four problem-specific callables (objective, move enumeration, move execution, tenure rule) from the memory logic. Use it to prototype; the two worked applications then show why production implementations specialize the scan and the delta bookkeeping for speed.

import numpy as np
from dataclasses import dataclass
from collections.abc import Callable, Hashable, Iterable
from typing import Any

@dataclass(frozen=True)
class Move:
    """One candidate move: payload to execute, objective change, tabu attributes."""
    payload: Any
    delta: float
    check_key: Hashable    # move is tabu while this attribute is active
    record_key: Hashable   # attribute activated when the move is executed

def tabu_search(
    initial: Any,
    objective: Callable[[Any], float],
    enumerate_moves: Callable[[Any], Iterable[Move]],
    execute: Callable[[Any, Any], Any],
    tenure: Callable[[np.random.Generator], int],
    n_iters: int,
    seed: int = 0,
) -> tuple[Any, float, list[float]]:
    """Attribute-based best-improvement tabu search (minimization).

    Returns (best solution, best value, best-so-far trace per iteration)."""
    rng = np.random.default_rng(seed)
    current, cur_val = initial, objective(initial)
    best, best_val = current, cur_val
    tabu_until: dict[Hashable, int] = {}
    trace = [best_val]
    for t in range(1, n_iters + 1):
        chosen, chosen_val = None, np.inf
        fallback, fallback_exp = None, np.inf   # least-tabu default aspiration
        for mv in enumerate_moves(current):
            new_val = cur_val + mv.delta
            if tabu_until.get(mv.check_key, 0) >= t and new_val >= best_val:
                if tabu_until[mv.check_key] < fallback_exp:
                    fallback, fallback_exp = mv, tabu_until[mv.check_key]
                continue
            if new_val < chosen_val:
                chosen, chosen_val = mv, new_val
        if chosen is None:                      # every move tabu, none aspired
            if fallback is None:
                break                           # neighborhood empty: stop
            chosen, chosen_val = fallback, cur_val + fallback.delta
        current = execute(current, chosen.payload)
        cur_val = chosen_val
        tabu_until[chosen.record_key] = t + tenure(rng)
        if cur_val < best_val:
            best, best_val = current, cur_val
        trace.append(best_val)
    return best, best_val, trace

# Demo: QUBO  min x^T Q x  with bit-flip moves and O(n) vectorized deltas.
def qubo_flip_deltas(Q: np.ndarray, x: np.ndarray) -> np.ndarray:
    """Objective change of flipping each bit of x (Q symmetric)."""
    s, dg = Q @ x, np.diag(Q)
    return (1 - 2 * x) * (dg + 2 * (s - dg * x))

rng = np.random.default_rng(7)
n = 14
Q = rng.integers(-10, 11, (n, n))
Q = (Q + Q.T) // 2

def moves(x: np.ndarray) -> list[Move]:
    """All bit flips; attribute = flipped index (element attribute)."""
    d = qubo_flip_deltas(Q, x)
    return [Move(payload=i, delta=float(d[i]), check_key=i, record_key=i)
            for i in range(n)]

def do_flip(x: np.ndarray, i: int) -> np.ndarray:
    """Return a copy of x with bit i flipped."""
    y = x.copy()
    y[i] = 1 - y[i]
    return y

x0 = rng.integers(0, 2, n)
best, best_val, trace = tabu_search(
    x0, lambda x: float(x @ Q @ x), moves, do_flip,
    tenure=lambda g: int(g.integers(5, 9)), n_iters=200, seed=1)
bits = ((np.arange(2 ** n)[:, None] >> np.arange(n)) & 1).astype(np.int64)
exact = int(np.einsum("bi,ij,bj->b", bits, Q, bits).min())
print(best_val, exact)
# Expected: tabu search reaches -143.0, equal to the brute-force optimum -143

Parameter guidance

Operator-level details (swap vs insertion neighborhoods, delta derivations per move type, perturbation design) belong to local-search-and-neighborhoods and fitness-evaluation-and-caching — this skill assumes a working neighborhood and adds memory on top.

Worked Application 1: Robust Tabu Search for the QAP

The quadratic assignment problem: given flow matrix $F \in \mathbb{Z}^{n \times n}$ and distance matrix $D \in \mathbb{Z}^{n \times n}$, find a permutation $\pi$ (unit $u$ placed at location $\pi(u)$) minimizing

$$ C(\pi) = \sum_{u=1}^{n} \sum_{v=1}^{n} F_{uv} , D_{\pi(u)\pi(v)}. $$

Exact methods stall around $n \approx 30$; Taillard's robust tabu search (Ro-TS) remains the reference heuristic and the standard baseline on QAPLIB (see quadratic-assignment-problem for the problem-side view). The swap move exchanges the locations of units $r$ and $s$. With $p = \pi(r)$, $q = \pi(s)$, the cost change is computable in O(n):

$$ \Delta(\pi, r, s) = F_{rr}(D_{qq}{-}D_{pp}) + F_{rs}(D_{qp}{-}D_{pq}) + F_{sr}(D_{pq}{-}D_{qp}) + F_{ss}(D_{pp}{-}D_{qq}) + \sum_{k \neq r,s} \big[ F_{kr}(D_{\pi(k)q}{-}D_{\pi(k)p}) + F_{ks}(D_{\pi(k)p}{-}D_{\pi(k)q}) + F_{rk}(D_{q\pi(k)}{-}D_{p\pi(k)}) + F_{sk}(D_{p\pi(k)}{-}D_{q\pi(k)}) \big] $$

Ro-TS maintains the full $\Delta$ matrix. After executing swap $(r, s)$, every pair ${u, v}$ disjoint from ${r, s}$ updates in O(1) (Taillard 1991), with $\pi$ taken before the swap:

$$ \Delta'(u,v) = \Delta(u,v) + (F_{ru}{-}F_{rv}{+}F_{sv}{-}F_{su})(D_{\pi(r)\pi(u)}{-}D_{\pi(r)\pi(v)}{+}D_{\pi(s)\pi(v)}{-}D_{\pi(s)\pi(u)}) + (F_{ur}{-}F_{vr}{+}F_{vs}{-}F_{us})(D_{\pi(u)\pi(r)}{-}D_{\pi(v)\pi(r)}{+}D_{\pi(v)\pi(s)}{-}D_{\pi(u)\pi(s)}) $$

Only the 2n−3 pairs touching $r$ or $s$ need the O(n) formula again, so one iteration costs O(n²) total — scan, choice, and maintenance included — versus O(n³) for naive re-evaluation. Both formulas below handle general (asymmetric) matrices and are verified against brute-force recomputation.

import numpy as np

def qap_cost(F: np.ndarray, D: np.ndarray, perm: np.ndarray) -> int:
    """Full objective: sum_{u,v} F[u,v] * D[perm[u], perm[v]]."""
    return int(np.sum(F * D[np.ix_(perm, perm)]))

def swap_delta(F: np.ndarray, D: np.ndarray, perm: np.ndarray,
               r: int, s: int) -> int:
    """O(n) cost change of swapping perm[r] and perm[s] (general matrices)."""
    pr, ps = perm[r], perm[s]
    d = (F[r, r] * (D[ps, ps] - D[pr, pr])
         + F[r, s] * (D[ps, pr] - D[pr, ps])
         + F[s, r] * (D[pr, ps] - D[ps, pr])
         + F[s, s] * (D[pr, pr] - D[ps, ps]))
    idx = np.arange(len(perm))
    k = np.flatnonzero((idx != r) & (idx != s))
    pk = perm[k]
    d += int(np.sum(F[k, r] * (D[pk, ps] - D[pk, pr])
                    + F[k, s] * (D[pk, pr] - D[pk, ps])
                    + F[r, k] * (D[ps, pk] - D[pr, pk])
                    + F[s, k] * (D[pr, pk] - D[ps, pk])))
    return int(d)

def full_delta_matrix(F: np.ndarray, D: np.ndarray,
                      perm: np.ndarray) -> np.ndarray:
    """All swap deltas; O(n^3) but computed once at the start of the search."""
    n = len(perm)
    delta = np.zeros((n, n), dtype=np.int64)
    for r in range(n - 1):
        for s in range(r + 1, n):
            delta[r, s] = swap_delta(F, D, perm, r, s)
            delta[s, r] = delta[r, s]
    return delta

The search itself. Tabu attribute = (unit, location) pair: after moving unit $r$ away from location $\pi(r)$, the assignment "$r$ back at $\pi(r)$" is tabu for a tenure drawn uniformly from $[0.9n, 1.1n]$ — redrawn for every recording, which is what makes the method robust to instance-specific tuning. A swap is tabu only if both units would re-enter recently-left locations. Two aspiration tests: improved-best, and the forced long-term test (an assignment absent for $5n^2$ iterations forces the move) — the latter is Ro-TS's built-in frequency-style diversification.

import numpy as np
# Continues the same module: uses qap_cost, swap_delta, full_delta_matrix
# from the previous block.

def robust_tabu_search(F: np.ndarray, D: np.ndarray, n_iters: int = 10_000,
                       seed: int = 0,
                       tenure: tuple[int, int] | None = None,
                       aspiration: int | None = None) -> tuple[np.ndarray, int]:
    """Taillard (1991) robust tabu search for the QAP; O(n^2) per iteration.

    Returns (best permutation, best cost)."""
    n = F.shape[0]
    rng = np.random.default_rng(seed)
    t_lo, t_hi = tenure if tenure else (int(0.9 * n), int(np.ceil(1.1 * n)))
    asp = aspiration if aspiration else 5 * n * n
    perm = rng.permutation(n)
    cost = qap_cost(F, D, perm)
    best_perm, best_cost = perm.copy(), cost
    delta = full_delta_matrix(F, D, perm)
    tabu_until = np.zeros((n, n), dtype=np.int64)      # [unit, location]
    iu, ju = np.triu_indices(n, 1)                     # all swaps, vectorized
    for t in range(1, n_iters + 1):
        d = delta[iu, ju]
        ti = tabu_until[iu, perm[ju]]   # unit iu entering location perm[ju]
        tj = tabu_until[ju, perm[iu]]   # unit ju entering location perm[iu]
        tabu = (ti >= t) & (tj >= t)    # tabu only if BOTH re-entries are recent
        aspired = (cost + d < best_cost) | (ti < t - asp) | (tj < t - asp)
        admissible = ~tabu | aspired
        if admissible.any():
            cand = np.flatnonzero(admissible)
            k = int(cand[np.argmin(d[cand])])
        else:
            k = int(np.argmin(d))       # all moves tabu: take the best anyway
        r, s = int(iu[k]), int(ju[k])
        cost += int(delta[r, s])
        a = F[r, :] - F[s, :]           # O(1)-per-pair update, pre-swap perm
        b = D[perm[r], perm] - D[perm[s], perm]
        c = F[:, r] - F[:, s]
        e = D[perm, perm[r]] - D[perm, perm[s]]
        delta += (np.subtract.outer(a, a) * np.subtract.outer(b, b)
                  + np.subtract.outer(c, c) * np.subtract.outer(e, e))
        tabu_until[r, perm[r]] = t + rng.integers(t_lo, t_hi + 1)
        tabu_until[s, perm[s]] = t + rng.integers(t_lo, t_hi + 1)
        perm[r], perm[s] = perm[s], perm[r]
        for u in range(n):              # exact O(n) refresh of pairs touching r, s
            for v in (r, s):
                if u != v:
                    dv = swap_delta(F, D, perm, u, v)
                    delta[u, v] = dv
                    delta[v, u] = dv
        if cost < best_cost:
            best_cost, best_perm = cost, perm.copy()
    return best_perm, best_cost

# Tiny synthetic instance, checked against brute force.
from itertools import permutations
rng = np.random.default_rng(42)
n = 8
F = rng.integers(0, 10, (n, n)); np.fill_diagonal(F, 0)
D = rng.integers(1, 10, (n, n)); np.fill_diagonal(D, 0)
opt = min(qap_cost(F, D, np.array(p)) for p in permutations(range(n)))
best_perm, best_cost = robust_tabu_search(F, D, n_iters=2000, seed=1)
print(opt, best_cost, qap_cost(F, D, best_perm))
# Expected: 1155 1155 1155 — Ro-TS matches the brute-force optimum, and the
# independently recomputed cost confirms the incremental bookkeeping.

Implementation notes that matter in practice:

• The incremental cost must be advanced with delta[r, s] before the outer-product update touches the matrix, and the update term must use the pre-swap permutation. Both orderings are easy to get wrong; the final qap_cost recomputation in the demo is the regression check that catches them.
• The "both re-entries recent" tabu test (logical AND) is deliberately weaker than OR. With OR, after a few iterations almost every swap is tabu and the search runs on aspiration overrides alone.
• For symmetric instances with zero diagonals the delta formula collapses to $2\sum_{k\neq r,s}(F_{sk}-F_{rk})(D_{\pi(r)\pi(k)}-D_{\pi(s)\pi(k)})$; keep the general form unless profiling says otherwise.

Worked Application 2: Job-Shop Scheduling with a Critical-Path Neighborhood

The job-shop problem (JSP): $n$ jobs, each a fixed sequence of operations over $m$ machines, one operation per machine per job; minimize makespan $C_{\max}$. A solution is an operation order for each machine. In the disjunctive-graph view, fixing all machine orders turns the schedule into a longest-path computation: each operation's start time is the longest path from the source, and $C_{\max}$ is the longest path overall. A critical path is a chain of operations realizing $C_{\max}$, linked by job-precedence or machine-precedence arcs.

Two classic results make the neighborhood (van Laarhoven, Aarts & Lenstra 1992, "Job shop scheduling by simulated annealing"):

1. Swapping two adjacent critical operations on the same machine never creates a cycle — every such move is feasible without checking.
2. Swapping adjacent operations that are not both critical can never reduce the makespan — so restricting moves to critical machine-arcs loses nothing. This is a candidate list strategy that is exact, not approximate.

The tabu attribute is the destroyed disjunctive arc: after swapping $u$ before $v$ into $v$ before $u$ on machine $M$, re-creating arc $u \to v$ on $M$ is tabu. This is the attribute family behind the best JSP tabu searches (Taillard 1994, "Parallel taboo search techniques for the job shop scheduling problem"; Nowicki & Smutnicki 1996, "A fast taboo search algorithm for the job shop problem"). Graph machinery first:

import numpy as np

def makespan(machines: np.ndarray, ptimes: np.ndarray,
             seq: list[list[int]]) -> tuple[int, np.ndarray, np.ndarray,
                                            np.ndarray, np.ndarray] | None:
    """Longest-path schedule for fixed machine sequences.

    machines[j, k] / ptimes[j, k]: machine and duration of job j's k-th op.
    seq[m]: job order on machine m. Operation id = j * n_mach + k.
    Returns (cmax, start, completion, job_pred, mach_pred); None if the
    sequences induce a cycle (only possible for externally supplied seq)."""
    n_jobs, n_mach = machines.shape
    n_ops = n_jobs * n_mach
    op_of = np.zeros((n_jobs, n_mach), dtype=np.int64)   # (job, machine) -> op
    for j in range(n_jobs):
        for k in range(n_mach):
            op_of[j, machines[j, k]] = j * n_mach + k
    p = ptimes.ravel()
    job_pred = np.full(n_ops, -1, dtype=np.int64)
    for j in range(n_jobs):
        job_pred[j * n_mach + 1: (j + 1) * n_mach] = \
            np.arange(j * n_mach, (j + 1) * n_mach - 1)
    mach_pred = np.full(n_ops, -1, dtype=np.int64)
    for m in range(n_mach):
        ops = op_of[seq[m], m]
        mach_pred[ops[1:]] = ops[:-1]
    succ_job = np.full(n_ops, -1, dtype=np.int64)
    succ_mach = np.full(n_ops, -1, dtype=np.int64)
    for o in range(n_ops):
        if job_pred[o] >= 0:
            succ_job[job_pred[o]] = o
        if mach_pred[o] >= 0:
            succ_mach[mach_pred[o]] = o
    indeg = (job_pred >= 0).astype(np.int64) + (mach_pred >= 0).astype(np.int64)
    start = np.zeros(n_ops, dtype=np.int64)
    comp = np.zeros(n_ops, dtype=np.int64)
    stack, seen = list(np.flatnonzero(indeg == 0)), 0
    while stack:                       # topological longest-path pass
        o = stack.pop()
        seen += 1
        st = 0
        if job_pred[o] >= 0:
            st = max(st, comp[job_pred[o]])
        if mach_pred[o] >= 0:
            st = max(st, comp[mach_pred[o]])
        start[o], comp[o] = st, st + p[o]
        for nx in (succ_job[o], succ_mach[o]):
            if nx >= 0:
                indeg[nx] -= 1
                if indeg[nx] == 0:
                    stack.append(nx)
    if seen < n_ops:
        return None
    return int(comp.max()), start, comp, job_pred, mach_pred

def critical_path(machines: np.ndarray, ptimes: np.ndarray,
                  seq: list[list[int]]) -> tuple[int, list[int], np.ndarray]:
    """One critical path (source to sink) and the machine-predecessor array."""
    res = makespan(machines, ptimes, seq)
    if res is None:
        raise ValueError("machine sequences induce a cycle")
    cmax, start, comp, job_pred, mach_pred = res
    o = int(np.argmax(comp))
    path = [o]
    while True:
        o = path[-1]
        if mach_pred[o] >= 0 and comp[mach_pred[o]] == start[o]:
            path.append(int(mach_pred[o]))
        elif job_pred[o] >= 0 and comp[job_pred[o]] == start[o]:
            path.append(int(job_pred[o]))
        else:
            break
    path.reverse()
    return cmax, path, mach_pred

def critical_machine_pairs(machines: np.ndarray, path: list[int],
                           mach_pred: np.ndarray) -> list[tuple[int, int, int]]:
    """Adjacent same-machine critical pairs as (machine, first_job, second_job)."""
    n_mach = machines.shape[1]
    pairs = []
    for a, b in zip(path[:-1], path[1:]):
        if mach_pred[b] == a:
            ja, jb = a // n_mach, b // n_mach
            pairs.append((int(machines[ja, a % n_mach]), int(ja), int(jb)))
    return pairs

The tabu search. Initial solution: order jobs on every machine by the operation's position inside its job (ties by job id) — this "level order" is provably acyclic and needs no scheduler. Tenure is randomized on a small range; the JSP arc attribute is much stronger than the QAP assignment attribute, so tenures of 8–14 work where QAP needs ~$n$. Each candidate move is evaluated by a full O(|ops|) recomputation, which is honest and simple; Taillard (1994) shows how head/tail values give an O(1) makespan lower bound per move when the scan itself becomes the bottleneck — see job-shop-scheduling for that refinement and for the Nowicki–Smutnicki block neighborhood.

import numpy as np
# Continues the same module: uses makespan, critical_path,
# critical_machine_pairs from the previous block.

def initial_sequences(machines: np.ndarray) -> list[list[int]]:
    """Acyclic start: order jobs on each machine by op position in the job."""
    n_jobs, n_mach = machines.shape
    pos = np.zeros((n_jobs, n_mach), dtype=np.int64)
    for j in range(n_jobs):
        for k in range(n_mach):
            pos[j, machines[j, k]] = k
    return [list(np.lexsort((np.arange(n_jobs), pos[:, m])))
            for m in range(n_mach)]

def apply_swap(seq: list[list[int]], m: int, ja: int, jb: int) -> list[list[int]]:
    """Copy of seq with jobs ja, jb (adjacent on machine m, ja first) swapped."""
    new = [list(s) for s in seq]
    i = new[m].index(ja)
    new[m][i], new[m][i + 1] = jb, ja
    return new

def jsp_tabu_search(machines: np.ndarray, ptimes: np.ndarray,
                    n_iters: int = 500, tenure: tuple[int, int] = (8, 14),
                    seed: int = 0) -> tuple[list[list[int]], int]:
    """Critical-path tabu search for the JSP (N1 neighborhood).

    Tabu attribute: the destroyed machine arc (machine, first_job, second_job).
    Returns (best machine sequences, best makespan)."""
    rng = np.random.default_rng(seed)
    seq = initial_sequences(machines)
    cmax, path, mach_pred = critical_path(machines, ptimes, seq)
    best_seq, best_cmax = [list(s) for s in seq], cmax
    tabu: dict[tuple[int, int, int], int] = {}
    for t in range(1, n_iters + 1):
        pairs = critical_machine_pairs(machines, path, mach_pred)
        if not pairs:
            break                                  # cmax = max job length: optimal
        chosen, chosen_val = None, None            # (move, seq, schedule result)
        fallback, fallback_exp = None, None
        for m, ja, jb in pairs:
            cand_seq = apply_swap(seq, m, ja, jb)  # feasible by construction
            res = critical_path(machines, ptimes, cand_seq)
            val = res[0]
            if tabu.get((m, jb, ja), 0) >= t and val >= best_cmax:
                exp = tabu[(m, jb, ja)]            # tabu, no aspiration
                if fallback_exp is None or exp < fallback_exp:
                    fallback, fallback_exp = ((m, ja, jb), cand_seq, res), exp
                continue
            if chosen_val is None or val < chosen_val:
                chosen, chosen_val = ((m, ja, jb), cand_seq, res), val
        if chosen is None:                         # all tabu: least-tabu fallback
            chosen = fallback
        (m, ja, jb), seq, (cmax, path, mach_pred) = chosen
        tabu[(m, ja, jb)] = t + int(rng.integers(tenure[0], tenure[1] + 1))
        if cmax < best_cmax:
            best_cmax, best_seq = cmax, [list(s) for s in seq]
    return best_seq, best_cmax

# ft06 (Fisher & Thompson 1963), 6 jobs x 6 machines, known optimum 55.
machines = np.array([[2, 0, 1, 3, 5, 4], [1, 2, 4, 5, 0, 3],
                     [2, 3, 5, 0, 1, 4], [1, 0, 2, 3, 4, 5],
                     [2, 1, 4, 5, 0, 3], [1, 3, 5, 0, 4, 2]])
ptimes = np.array([[1, 3, 6, 7, 3, 6], [8, 5, 10, 10, 10, 4],
                   [5, 4, 8, 9, 1, 7], [5, 5, 5, 3, 8, 9],
                   [9, 3, 5, 4, 3, 1], [3, 3, 9, 10, 4, 1]])
start_cmax = critical_path(machines, ptimes, initial_sequences(machines))[0]
best_seq, best_cmax = jsp_tabu_search(machines, ptimes, n_iters=500, seed=0)
print(start_cmax, best_cmax)
# Expected: 60 55 — from the level-order start (makespan 60) the critical-path
# tabu search reaches 55, the proven optimum of ft06, within 500 iterations.

Two observations worth internalizing. First, the neighborhood here is the candidate list: instead of all O(n·m) adjacent swaps, only critical-arc swaps are scanned (typically 3–10 on ft06-sized instances), and result 2 above guarantees this filter discards no improving move. Second, the uphill steps are essential: from many schedules every critical swap increases the makespan, and a memoryless descent stops exactly there; the tabu list is what lets the search traverse those ridges without immediately undoing them.

Advanced Techniques

Reactive tenure

Battiti & Tecchiolli (1994) replace tenure tuning with feedback: hash every visited solution; if a solution repeats within a window, the search is cycling — increase tenure multiplicatively. If no repetition occurs for a long stretch, decay tenure so the search is no more restricted than necessary. Canonical hashing of solutions is covered in fitness-evaluation-and-caching.

import numpy as np

class ReactiveTenure:
    """Battiti & Tecchiolli (1994): adapt tenure from solution repetitions."""

    def __init__(self, start: int, cap: int, window: int = 50,
                 up: float = 1.1, down: float = 0.9) -> None:
        self.tenure, self.cap, self.window = start, cap, window
        self.up, self.down = up, down
        self.seen: dict[int, int] = {}     # solution hash -> last iteration
        self.last_change = 0

    def update(self, solution_hash: int, t: int) -> int:
        """Call once per iteration with the current solution's hash."""
        last = self.seen.get(solution_hash)
        self.seen[solution_hash] = t
        if last is not None and t - last < self.window:        # cycling
            self.tenure = min(self.cap, max(self.tenure + 1,
                                            int(self.tenure * self.up)))
            self.last_change = t
        elif t - self.last_change > self.window:               # quiet: relax
            self.tenure = max(1, int(self.tenure * self.down))
            self.last_change = t
        return self.tenure

rt = ReactiveTenure(start=8, cap=40)
perm = np.array([3, 1, 0, 2])
h = hash(tuple(perm))
print(rt.update(h, t=1), rt.update(h, t=5))
# Expected: 8 9 — the repeat at t=5 (within the window) raises the tenure.

Frequency-based long-term memory and diversification

Recency memory forgets after $\tau$ iterations; frequency memory never forgets. Track transition counts (how often each move attribute fired) and residence counts (how long each attribute stayed in the solution). When the incumbent stalls, enter a diversification phase: add a penalty $w \cdot \text{freq}(m)/t$ to each move's delta so the scan prefers rarely used moves, or force the least-used assignment outright (the Ro-TS forced aspiration is exactly this, embedded in the aspiration test). Scale $w$ to the median |delta| of recently accepted moves so the penalty competes with, but does not drown, the objective.

import numpy as np

class FrequencyMemory:
    """Long-term transition/residence counts for a permutation problem."""

    def __init__(self, n: int) -> None:
        self.transition = np.zeros((n, n), dtype=np.int64)   # swap (i, j) count
        self.residence = np.zeros((n, n), dtype=np.int64)    # iters of i at pos
        self.iters = 0

    def record(self, perm: np.ndarray, r: int, s: int) -> None:
        """Call once per iteration, after executing swap (r, s)."""
        self.transition[r, s] += 1
        self.transition[s, r] += 1
        self.residence[np.arange(len(perm)), perm] += 1
        self.iters += 1

    def penalties(self, weight: float) -> np.ndarray:
        """Penalty matrix to ADD to the delta matrix while diversifying."""
        if self.iters == 0:
            return np.zeros_like(self.transition, dtype=np.float64)
        return weight * self.transition / self.iters

mem = FrequencyMemory(5)
perm = np.array([2, 0, 1, 4, 3])
mem.record(perm, 0, 1)
mem.record(perm, 0, 1)
mem.record(perm, 2, 3)
print(round(mem.penalties(10.0)[0, 1], 2), round(mem.penalties(10.0)[2, 3], 2))
# Expected: 6.67 3.33 — the over-used swap (0, 1) is penalized twice as hard.

Intensification by restarting from elite solutions

Quality memory: keep a bounded pool of the best distinct solutions seen. When diversification has not produced a new incumbent for a long phase, restart the search from a random elite with a cleared (or shortened) tabu list and a fresh tenure draw — the second visit explores different uphill exits because the randomized tenure changes which reversals are forbidden. Nowicki & Smutnicki's TSAB owes much of its strength to exactly this back-jump-to-elite mechanism. A stronger variant fixes the attributes shared by all elites (common assignments, common arcs) and searches only the residual problem; path relinking between elites is the systematic version of that idea. Restart policy comparisons and diversity measurement belong to diversity-and-population-management.

import numpy as np

class ElitePool:
    """Bounded pool of best distinct solutions for intensification restarts."""

    def __init__(self, capacity: int) -> None:
        self.capacity = capacity
        self.pool: list[tuple[float, tuple[int, ...]]] = []

    def offer(self, value: float, solution: np.ndarray) -> None:
        """Insert unless duplicate; keep the `capacity` best by value."""
        key = tuple(int(v) for v in solution)
        if any(k == key for _, k in self.pool):
            return
        self.pool.append((value, key))
        self.pool.sort(key=lambda e: e[0])
        del self.pool[self.capacity:]

    def restart_from(self, rng: np.random.Generator) -> np.ndarray:
        """Pick an elite uniformly at random as the new current solution."""
        _, sol = self.pool[rng.integers(len(self.pool))]
        return np.array(sol)

pool = ElitePool(capacity=3)
for v, p in [(10, [0, 1, 2]), (8, [1, 0, 2]), (12, [2, 1, 0]),
             (8, [1, 0, 2]), (7, [2, 0, 1])]:
    pool.offer(float(v), np.array(p))
print([v for v, _ in pool.pool])
# Expected: [7.0, 8.0, 10.0] — duplicate rejected, worst (12) evicted.

Candidate list strategies

When the full neighborhood scan is too slow, scan a subset — but choose the subset with structure, not blindly (Glover & Laguna 1997, ch. 3):

• Exact filters: problem theory proves some moves cannot improve — the JSP critical-arc restriction above is the model case. Always look for one first.
• Elite candidate list: scan the full neighborhood occasionally, keep the best $k$ moves; in between, scan only that list, refreshing it when its best move degrades below a threshold.
• Successive filtering: cheap test first (e.g., lower bound on delta), exact evaluation only for survivors.
• Random sampling: evaluate a random fraction per iteration; simplest, weakest, but unbiased — combine with first-improvement-style early exit.
• Neighborhood decomposition: partition moves (by region, machine, vehicle) and rotate through partitions; pairs naturally with frequency memory to decide which partition deserves attention.

Strategic oscillation

For tightly constrained problems, let the search cross the feasibility boundar

…(truncated)