Biased Random Key Genetic Algorithm

Biased Random-Key Genetic Algorithm (BRKGA)

You are an expert in biased random-key genetic algorithms for combinatorial optimization. This skill covers the random-key encoding, the elite/mutant population partition, biased uniform crossover, decoder design as the single problem-specific component, and a reusable numpy framework with two complete worked applications (permutation flow shop with a sort decoder, set covering with a threshold decoder). Use the framework below to assess whether BRKGA fits the problem, build a correct and efficient implementation, and diagnose convergence problems.

Initial Assessment

Before writing any BRKGA code, establish the following:

• Solution structure. What object does a solution decode to: a permutation, a subset, an assignment, a schedule, or a combination? This single fact determines the decoder family (sort, threshold, greedy-priority, multi-segment) and therefore the chromosome length.
• Chromosome length n. Count the decisions one key must drive. For sequencing, n = number of jobs/items. For selection, n = number of candidate elements. Multi-segment chromosomes concatenate one key block per decision layer.
• Decode cost. Decoding dominates BRKGA runtime; the genetic operators are O(p·n) memory copies and never the bottleneck. Estimate the cost of one decode and multiply by pop_size × generations. If a single decode takes more than a few milliseconds, plan for vectorized batch decoding, caching of elite fitness, or parallel decoding from the start.
• Constraint placement. Decide, per constraint family, whether the decoder absorbs it (construct only feasible solutions), repairs it (fix violations after a cheap construction), or penalizes it (add a violation term to fitness). BRKGA's main selling point is that the decoder can guarantee feasibility, so prefer absorb/repair over penalties.
• Objective. BRKGA as described here minimizes a single scalar. Multi-objective variants exist (NSGA-II-style sorting on top of the BRKGA population) but need extra machinery.
• Time budget and stopping rule. Fixed generation count, wall-clock limit, or stall limit (no improvement for k generations)? This drives population size: with a tight budget, prefer a smaller population and more generations.
• Baseline and warm start. Is there a constructive heuristic (greedy, NEH, LPT) whose output can be encoded as keys and injected into the initial population? Always benchmark BRKGA against that baseline; a metaheuristic that loses to greedy is misconfigured.
• Library vs. from scratch. For production use, the maintained brkga_mp_ipr packages (Python/C++/Julia) implement multi-parent BRKGA with implicit path relinking. Implement from scratch (as below) when you need full control over the decoder/evaluation loop, vectorization across the population, or research instrumentation.
• Exact alternative. If instances are small enough for a MIP solver to close the gap in the available time, use the exact model and keep BRKGA for the large instances or as a warm-start provider.
• Reproducibility. Fix seeds for instance generation and for the algorithm separately. Plan multiple independent runs (different algorithm seeds) if you will report statistics.

Algorithm Anatomy

BRKGA (Gonçalves & Resende, 2011, "Biased random-key genetic algorithms for combinatorial optimization") combines three ideas on top of Bean's random-key GA (Bean, 1994, "Genetic algorithms and random keys for sequencing and optimization"):

1. Random-key encoding. A chromosome is a vector $x \in [0,1)^n$. A problem-specific decoder $D$ maps any such vector to a feasible solution, and fitness is

$$ \text{fit}(x) = f(D(x)), \qquad D : [0,1)^n \to S, $$

where $S$ is the feasible set. Because every point of the hypercube decodes to something feasible, all genetic operators act on $[0,1)^n$ and never require problem-specific repair. The decoder is the only part you rewrite per problem.

2. Elite/mutant partition. Each generation, sort the population of size $p$ by fitness and split it:

• Elite set $E$: the best $n_e = \lceil p_e \cdot p \rceil$ chromosomes, copied unchanged to the next generation (strong elitism; best fitness is monotone).
• Mutants: $n_m = \lceil p_m \cdot p \rceil$ fresh uniform random vectors. BRKGA has no per-gene mutation operator; mutants play that diversification role, like immigrants.
• Offspring: the remaining $p - n_e - n_m$ slots, produced by biased crossover.

3. Biased uniform crossover. Each offspring takes one parent uniformly from the elite set and one from the non-elite set. Gene $i$ comes from the elite parent with probability $\rho_e > 0.5$:

$$ c_i = \begin{cases} a_i & \text{with probability } \rho_e \quad (a \in E) \ b_i & \text{otherwise} \quad (b \notin E) \end{cases} $$

This is parameterized uniform crossover (Spears & De Jong, 1991, "On the virtues of parameterized uniform crossover") with the bias always pointing at the better parent. With $\rho_e = 0.7$ an offspring inherits 70% of its genes from an elite chromosome in expectation, so the search drifts toward elite regions while non-elite parents and mutants keep injecting variation. The requirement $p_e < 0.5$ keeps the elite pool smaller than the non-elite pool, which maintains selection pressure.

Decoder taxonomy. Choosing the decoder is the real design decision; see decoder-based-representations for the full catalog. The common families:

When to use BRKGA — and when not.

• Use BRKGA when feasibility is hard to preserve under direct crossover/mutation, when a priority-driven constructive heuristic already exists (wrap it as a greedy decoder), or when you want one tested evolutionary engine reused across many problems with only the decoder changing.
• Prefer a permutation GA with OX/PMX operators (see genetic-algorithms) when adjacency carries the fitness signal (pure TSP): sort decoders preserve relative order under crossover, not adjacency, so BRKGA needs a local-search hybrid to be competitive there.
• Prefer trajectory methods (tabu search, ILS) when one decode is expensive and a fast incremental (delta) evaluation of small moves exists — population methods cannot exploit delta evaluation as directly.

Complexity per generation. Sorting is $O(p \log p)$; crossover and mutant generation are $O(p n)$ vectorized array operations; decoding costs $(p - n_e) \cdot C_{dec}$ because elite fitness is cached. With a sort decoder, $C_{dec} = O(n \log n)$ plus the objective evaluation. The decoder is almost always the bottleneck — vectorize it across the population whenever the construction logic allows.

Landscape note. The decoder is many-to-one: uncountably many key vectors map to the same solution (for a sort decoder, only the ordering of keys matters). This redundancy creates plateaus in key space. It is harmless to the dynamics but means diversity must be measured on decoded solutions or fitness values, never on raw key distances.

Core Framework Implementation

The framework is fully problem-independent. Pseudocode first:

BRKGA(decode, p, p_e, p_m, rho_e, G):
    population <- p random vectors drawn uniformly from [0,1)^n
    fitness    <- decode(population)                      # batch evaluation
    repeat G generations:
        sort population by fitness (ascending; minimization)
        ELITE  <- best ceil(p_e * p) chromosomes          # copied unchanged
        MUTANT <- ceil(p_m * p) fresh uniform random vectors
        for each of the p - |ELITE| - |MUTANT| offspring slots:
            a <- chromosome drawn uniformly from ELITE
            b <- chromosome drawn uniformly from NON-ELITE
            for each gene i in 1..n:
                child[i] <- a[i] with probability rho_e, else b[i]
        population <- ELITE  ∪  offspring  ∪  MUTANT
        fitness    <- cached elite fitness  ∪  decode(offspring ∪ MUTANT)
    return best chromosome, decode(best)

Save the following as brkga_core.py; both worked examples import from it. The decode callable receives the whole key matrix at once so problem code can vectorize across the population, and elite fitness values are never recomputed.

"""brkga_core.py — generic BRKGA framework. Only the decoder is problem-specific."""
from dataclasses import dataclass
from typing import Callable

import numpy as np


@dataclass(frozen=True)
class BrkgaParams:
    """Configuration of one BRKGA run (minimization)."""
    n_genes: int
    pop_size: int = 100
    elite_frac: float = 0.20
    mutant_frac: float = 0.15
    rho_e: float = 0.70
    n_generations: int = 200
    seed: int = 0


def brkga(
    decode: Callable[[np.ndarray], np.ndarray],
    params: BrkgaParams,
    init_keys: np.ndarray | None = None,
) -> tuple[np.ndarray, float, list[float]]:
    """Run BRKGA; return (best chromosome, best fitness, best-per-generation history).

    `decode` maps a key matrix of shape (P, n_genes) to a fitness vector of
    shape (P,); lower is better. Rows of `init_keys` (if given) are injected
    into the initial population, e.g. to warm-start from a known solution.
    """
    rng = np.random.default_rng(params.seed)
    p, n = params.pop_size, params.n_genes
    n_elite = max(1, int(round(params.elite_frac * p)))
    n_mutant = max(1, int(round(params.mutant_frac * p)))
    n_offspring = p - n_elite - n_mutant

    pop = rng.random((p, n))
    if init_keys is not None:
        pop[: len(init_keys)] = init_keys
    fit = np.asarray(decode(pop), dtype=float)

    history: list[float] = []
    for _ in range(params.n_generations):
        order = np.argsort(fit)
        pop, fit = pop[order], fit[order]
        history.append(float(fit[0]))

        elite_parent = rng.integers(0, n_elite, size=n_offspring)
        other_parent = rng.integers(n_elite, p, size=n_offspring)
        from_elite = rng.random((n_offspring, n)) < params.rho_e
        offspring = np.where(from_elite, pop[elite_parent], pop[other_parent])
        mutants = rng.random((n_mutant, n))

        new_part = np.vstack([offspring, mutants])
        pop = np.vstack([pop[:n_elite], new_part])
        fit = np.concatenate([fit[:n_elite], np.asarray(decode(new_part), dtype=float)])

    best = int(np.argmin(fit))
    history.append(float(fit[best]))
    return pop[best].copy(), float(fit[best]), history

Implementation notes:

• The biased crossover is one np.where over a Boolean mask — all offspring for the generation are produced in a single vectorized expression.
• Elite chromosomes carry their cached fitness forward; only offspring + mutants are decoded each generation. With p_e = 0.2, that alone saves 20% of all decode work.
• The framework never inspects the decoded solutions. Re-decode the returned best chromosome in problem code to materialize the actual solution object.
• The decoder must be a deterministic function of the keys. If it were stochastic, cached elite fitness would diverge from re-decoded values and the monotonicity guarantee (and its use as a correctness check) would be lost.
• Maximization problems: negate the objective inside decode; keep the framework minimizing.

Parameter guidance (recommended ranges follow Gonçalves & Resende, 2011):

Selection, replacement, and crossover are fixed by the method — do not graft tournament selection or per-gene mutation onto BRKGA; if you need those degrees of freedom, you want a classic GA (see genetic-algorithms) instead.

Worked Example 1: Permutation Flow Shop with a Sort Decoder

Problem. $n$ jobs are processed by $m$ machines in the same machine order; one permutation $\pi$ of the jobs is used on every machine (permutation flow shop, $F_m \mid prmu \mid C_{\max}$). With processing times $p_{ij}$ for machine $i$ and job $j$, completion times follow the recurrence

$$ C_{i,\pi(k)} = \max!\left(C_{i,\pi(k-1)},; C_{i-1,\pi(k)}\right) + p_{i,\pi(k)}, $$

and the objective is the makespan $C_{\max} = C_{m,\pi(n)}$. The problem is NP-hard for $m \ge 3$. See flow-shop-scheduling for exact models, NEH, and the full heuristic landscape; here it serves as the canonical sort-decoder application.

Encoding. One key per job; argsort(keys) is the job sequence. Crossover on keys can produce any permutation, and every key vector decodes to a valid sequence — no repair, ever.

The instance generator and a population-batched makespan evaluator (the recurrence runs once over jobs and machines, but for all chromosomes simultaneously):

"""flow_shop.py — instance generator and batched makespan evaluation."""
import numpy as np


def random_flow_shop(n_jobs: int, n_machines: int, seed: int) -> np.ndarray:
    """Taillard-style instance: integer times in [1, 99], shape (n_machines, n_jobs)."""
    rng = np.random.default_rng(seed)
    return rng.integers(1, 100, size=(n_machines, n_jobs)).astype(float)


def makespan_batch(perms: np.ndarray, proc: np.ndarray) -> np.ndarray:
    """Makespan of each job permutation. perms: (P, n) int; proc: (m, n) float."""
    n_pop, n_jobs = perms.shape
    n_mach = proc.shape[0]
    seq_times = proc[:, perms]                # (m, P, n): times in sequence order
    comp = np.zeros((n_pop, n_mach))          # completion time per (chromosome, machine)
    for j in range(n_jobs):
        stage = seq_times[:, :, j].T          # (P, m): j-th sequenced job, each machine
        comp[:, 0] += stage[:, 0]
        for i in range(1, n_mach):
            comp[:, i] = np.maximum(comp[:, i], comp[:, i - 1]) + stage[:, i]
    return comp[:, -1]

The loop nest is $O(nm)$ Python iterations, but each iteration is vectorized over the whole population, so evaluating 100 chromosomes costs barely more than evaluating one. The decoder and driver:

"""run_flow_shop.py — BRKGA with a sort decoder on the permutation flow shop."""
import numpy as np

from brkga_core import BrkgaParams, brkga
from flow_shop import makespan_batch, random_flow_shop


def make_flow_shop_decoder(proc: np.ndarray):
    """Sort decoder: argsort of the keys is the job sequence."""
    def decode(keys: np.ndarray) -> np.ndarray:
        perms = np.argsort(keys, axis=1)      # (P, n) permutations, vectorized
        return makespan_batch(perms, proc)
    return decode


def makespan_single(perm: np.ndarray, proc: np.ndarray) -> float:
    """Independent scalar validator for one permutation."""
    n_mach = proc.shape[0]
    comp = np.zeros(n_mach)
    for j in perm:
        comp[0] += proc[0, j]
        for i in range(1, n_mach):
            comp[i] = max(comp[i], comp[i - 1]) + proc[i, j]
    return float(comp[-1])


if __name__ == "__main__":
    proc = random_flow_shop(n_jobs=20, n_machines=5, seed=42)
    params = BrkgaParams(n_genes=20, pop_size=100, elite_frac=0.20,
                         mutant_frac=0.15, rho_e=0.70, n_generations=200, seed=1)
    best_keys, best_makespan, history = brkga(make_flow_shop_decoder(proc), params)
    best_perm = np.argsort(best_keys)

    assert makespan_single(best_perm, proc) == best_makespan
    print(f"makespan: {history[0]:.0f} (gen 0 best) -> {best_makespan:.0f} (final)")
    print("sequence:", best_perm.tolist())
    # Expected: makespan 1351 (gen 0 best) -> 1243 (final) for these exact seeds;
    # history is non-increasing because the elite set is copied unchanged.

Two things to verify on any sort-decoder application: (a) the validator (makespan_single) recomputes the objective independently of the batched evaluator — they must agree exactly; (b) the history is monotone non-increasing — if it is not, elite fitness caching or the sort is broken.

Worked Example 2: Set Covering with a Threshold Decoder

Problem. Given a 0/1 matrix $A \in {0,1}^{r \times c}$ and column costs $c_j > 0$, select a minimum-cost set of columns covering every row:

$$ \min \sum_{j} c_j y_j \quad \text{s.t.} \quad \sum_{j} a_{ij}, y_j \ge 1 ;; \forall i, \qquad y \in {0,1}^c. $$

Encoding. One key per column. The threshold decoder selects column $j$ when $x_j \ge 0.5$. Unlike the sort decoder, the raw selection can be infeasible (uncovered rows) or wasteful (redundant columns), so the decoder embeds two deterministic post-steps:

1. Greedy repair — while rows remain uncovered, add the unselected column minimizing cost per newly covered row (the classic Chvátal greedy ratio).
2. Redundancy pruning — scan selected columns in decreasing cost order and drop any column whose rows are all covered at least twice.

Because repair and pruning are deterministic, the decoder is still a well-defined function from keys to feasible solutions. Instance generator first:

"""set_cover.py (part 1) — random SCP instance generator."""
import numpy as np


def random_set_cover(
    n_rows: int, n_cols: int, density: float, seed: int
) -> tuple[np.ndarray, np.ndarray]:
    """Random SCP: boolean cover matrix A (n_rows x n_cols) and integer column costs.

    Patches the random matrix so every row is coverable and no column is empty,
    which guarantees feasibility of the instance.
    """
    rng = np.random.default_rng(seed)
    A = rng.random((n_rows, n_cols)) < density
    for i in np.flatnonzero(~A.any(axis=1)):      # every row coverable
        A[i, rng.integers(n_cols)] = True
    for j in np.flatnonzero(~A.any(axis=0)):      # no useless columns
        A[rng.integers(n_rows), j] = True
    cost = rng.integers(1, 100, size=n_cols).astype(float)
    return A, cost

Repair, pruning, and an independent feasibility validator:

"""set_cover.py (part 2) — repair, redundancy pruning, and validation."""
import numpy as np


def repair_and_prune(selected: np.ndarray, A: np.ndarray, cost: np.ndarray) -> np.ndarray:
    """Greedy repair (min cost per newly covered row), then drop redundant columns."""
    selected = selected.copy()
    covered = A[:, selected].any(axis=1) if selected.any() else np.zeros(A.shape[0], bool)
    while not covered.all():
        gain = A[~covered].sum(axis=0).astype(float)   # newly covered rows per column
        gain[selected] = 0.0
        ratio = np.where(gain > 0, cost / np.maximum(gain, 1e-12), np.inf)
        j = int(np.argmin(ratio))
        selected[j] = True
        covered |= A[:, j]
    count = A[:, selected].sum(axis=1)                 # coverage multiplicity per row
    for j in sorted(np.flatnonzero(selected), key=lambda col: -cost[col]):
        rows_j = A[:, j]
        if (count[rows_j] >= 2).all():
            selected[j] = False
            count[rows_j] -= 1
    return selected


def validate_cover(selected: np.ndarray, A: np.ndarray, cost: np.ndarray) -> float:
    """Independent check: raise if any row is uncovered, else return total cost."""
    if not A[:, selected].any(axis=1).all():
        raise ValueError("infeasible: some row is uncovered")
    return float(cost[selected].sum())

Decoder and driver. The thresholding is vectorized; repair is inherently sequential (each greedy pick changes the next pick), so it runs per chromosome — this is the expected shape for repair-based decoders, and the reason decode cost must be assessed up front:

"""run_set_cover.py — BRKGA with a threshold decoder on set covering."""
import numpy as np

from brkga_core import BrkgaParams, brkga
from set_cover import random_set_cover, repair_and_prune, validate_cover


def make_scp_decoder(A: np.ndarray, cost: np.ndarray):
    """Threshold decoder: key >= 0.5 selects a column; repair restores feasibility."""
    def decode(keys: np.ndarray) -> np.ndarray:
        raw = keys >= 0.5                       # vectorized over the population
        vals = np.empty(len(keys))
        for k in range(len(keys)):              # repair is sequential per chromosome
            sel = repair_and_prune(raw[k], A, cost)
            vals[k] = cost[sel].sum()
        return vals
    return decode


if __name__ == "__main__":
    A, cost = random_set_cover(n_rows=40, n_cols=120, density=0.06, seed=7)
    params = BrkgaParams(n_genes=120, pop_size=120, elite_frac=0.20,
                         mutant_frac=0.15, rho_e=0.70, n_generations=150, seed=1)
    best_keys, best_cost, history = brkga(make_scp_decoder(A, cost), params)
    best_sel = repair_and_prune(best_keys >= 0.5, A, cost)

    greedy = repair_and_prune(np.zeros(120, dtype=bool), A, cost)
    print(f"BRKGA cost {validate_cover(best_sel, A, cost):.0f} "
          f"({int(best_sel.sum())} columns), greedy {validate_cover(greedy, A, cost):.0f}")
    # Expected: BRKGA cost 361 (17 columns) vs greedy 392 for these exact seeds;
    # BRKGA's gen-0 best is 417 and improves monotonically to 361.

Note the built-in baseline: decoding the all-zeros selection runs pure greedy-plus-pruning. BRKGA must beat it; if it does not, the decoder, parameters, or budget are wrong. Also note the fitness reported by decode is the repaired solution's cost — the chromosome itself is not changed (a "Cartesian" decoder). The Lamarckian alternative that writes the repair back into the keys is covered under Advanced Techniques.

Advanced Techniques

Warm starts and solution injection

Any known solution can be encoded as keys and injected via init_keys. For a sort decoder, build keys whose argsort reproduces a given permutation (e.g., from NEH — Nawaz, Enscore & Ham, 1983 — for flow shop, or LPT for parallel-machine-scheduling dispatch decoders):

import numpy as np


def keys_from_permutation(perm: np.ndarray) -> np.ndarray:
    """Keys whose argsort reproduces `perm`; values spread uniformly over (0, 1)."""
    n = len(perm)
    keys = np.empty(n)
    keys[perm] = (np.arange(n) + 0.5) / n
    return keys


def jittered_copies(keys: np.ndarray, n_copies: int, sigma: float, seed: int) -> np.ndarray:
    """Stack of perturbed copies, clipped to [0, 1); avoids seeding identical clones."""
    rng = np.random.default_rng(seed)
    noise = rng.normal(0.0, sigma, size=(n_copies, len(keys)))
    return np.clip(keys + noise, 0.0, np.nextafter(1.0, 0.0))


perm = np.array([3, 0, 2, 1])
seeds = jittered_copies(keys_from_permutation(perm), n_copies=3, sigma=0.02, seed=0)
assert (np.argsort(keys_from_permutation(perm)) == perm).all()
# Expected: assertion passes; `seeds` holds 3 near-copies that mostly decode to perm.

Inject one exact copy plus a few jittered copies — never fill a large share of the population with one solution, or the elite set collapses to clones in the first generation. See warm-starts discussion in the Practical Challenges section below for the failure mode.

Multi-parent crossover and implicit path relinking (BRKGA-MP-IPR)

Andrade, Toso, Gonçalves & Resende (2021, "The Multi-Parent Biased Random-Key Genetic Algorithm with Implicit Path Relinking and its real-world applications") generalize the two-parent crossover: each offspring has $\pi_t$ parents ($\pi_e$ of them elite), and each gene is copied from a parent chosen with rank-based bias weights, e.g. $w_r \propto 1/r$ for parent rank $r$. Two parents with $\rho_e = w_1/(w_1+w_2)$ recovers classic BRKGA.

import numpy as np


def multi_parent_offspring(
    pop_sorted: np.ndarray,
    n_elite: int,
    n_offspring: int,
    total_parents: int,
    elite_parents: int,
    rng: np.random.Generator,
) -> np.ndarray:
    """BRKGA-MP gene transfer: each gene copied from one of `total_parents` parents,
    chosen with rank bias w_r = 1/r (Andrade et al., 2021). Assumes fitness-sorted pop."""
    p, n = pop_sorted.shape
    e_idx = rng.integers(0, n_elite, size=(n_offspring, elite_parents))
    o_idx = rng.integers(n_elite, p, size=(n_offspring, total_parents - elite_parents))
    parents = np.concatenate([e_idx, o_idx], axis=1)
    ranked = np.sort(parents, axis=1)          # lower population index = better rank
    weights = 1.0 / np.arange(1, total_parents + 1)
    choice = rng.choice(total_parents, size=(n_offspring, n), p=weights / weights.sum())
    src = ranked[np.arange(n_offspring)[:, None], choice]
    return pop_sorted[src, np.arange(n)]


rng = np.random.default_rng(0)
pop = rng.random((20, 6))
children = multi_parent_offspring(pop, n_elite=4, n_offspring=10,
                                  total_parents=5, elite_parents=2, rng=rng)
# Expected: children.shape == (10, 6); every gene equals some parent's gene.

Implicit path relinking (IPR) intensifies between two elite chromosomes without a problem-specific neighborhood: walk from one key vector toward the other by adopting blocks of keys (direct IPR) or block-induced permutation fragments (permutation IPR), decoding intermediate points and keeping the best. It is the standard answer when plain BRKGA stalls a few percent above best-known values; use the brkga_mp_ipr package rather than reimplementing IPR.

Shaking and restarts

BRKGA's elitism makes stagnation visible: the best value flatlines. Two standard escapes, in increasing strength — shaking re-randomizes a fraction of elite genes in place; a restart rebuilds the population from scratch, optionally re-injecting a shaken incumbent:

import dataclasses

import numpy as np

from brkga_core import BrkgaParams, brkga


def shake(elite_keys: np.ndarray, intensity: float, rng: np.random.Generator) -> np.ndarray:
    """Re-randomize a fraction `intensity` of elite genes (Andrade et al., 2021)."""
    mask = rng.random(elite_keys.shape) < intensity
    return np.where(mask, rng.random(elite_keys.shape), elite_keys)


def brkga_multistart(decode, params: BrkgaParams, n_starts: int,
                     shake_intensity: float = 0.2, n_injected: int = 2):
    """Sequential restarts; each restart re-seeds the population and injects
    shaken copies of the incumbent best chromosome."""
    rng = np.random.default_rng(params.seed)
    best_keys, best_val = None, np.inf
    for s in range(n_starts):
        run = dataclasses.replace(params, seed=params.seed + 1000 * (s + 1))
        init = None
        if best_keys is not None:
            init = shake(np.tile(best_keys, (n_injected, 1)), shake_intensity, rng)
        keys, val, _ = brkga(decode, run, init_keys=init)
        if val < best_val:
            best_keys, best_val = keys.copy(), val
    return best_keys, best_val
# Expected: returned best_val <= the single-run value for the same total budget
# on stagnation-prone instances; equal otherwise.

A practical trigger: restart when the best value has not improved for n_genes to 2 * n_genes generations. Spend a long time budget as several restarts rather than one long run.

Island model with elite migration

Run $K$ independent populations and periodically migrate elites along a ring. Islands restore diversity at the population level and parallelize trivially (one process per island):

import numpy as np


def migrate_ring(pops: list[np.ndarray], fits: list[np.ndarray], n_migrants: int) -> None:
    """Ring migration: best n_migrants of island k overwrite the worst of island k+1.
    Migrant fitness travels with the keys, so no re-decoding is needed."""
    snapshots = [
        (pop[np.argsort(fit)[:n_migrants]].copy(), np.sort(fit)[:n_migrants])
        for pop, fit in zip(pops, fits)
    ]
    for k in range(len(pops)):
        dst = (k + 1) % len(pops)
        worst = np.argsort(fits[dst])[-n_migrants:]
        pops[dst][worst], fits[dst][worst] = snapshots[k]


rng = np.random.default_rng(0)
pops = [rng.random((30, 8)) for _ in range(4)]
fits = [rng.random(30) for _ in range(4)]
migrate_ring(pops, fits, n_migrants=2)
# Expected: each island now contains the previous island's two best chromosomes.

Typical settings: 2–4 islands, migrate every 25–100 generations, 1–3 migrants. More frequent migration makes the islands behave like one population and loses the benefit.

Lamarckian decoders and online parameter control

A Lamarckian decoder writes its repair/local-search result back into the chromosome so the improvement is inherited. For a threshold decoder, reflect each key to the side of 0.5 matching the repaired solution — deterministic and information-preserving:

import numpy as np


def lamarckian_threshold_keys(keys: np.ndarray, selected: np.ndarray) -> np.ndarray:
    """Push each key to the side of 0.5 that matches the repaired selection."""
    flip = selected != (keys >= 0.5)
    new = keys.copy()
    new[flip] = 1.0 - keys[flip]
    return new


keys = np.array([0.10, 0.90, 0.60, 0.40])
selected = np.array([True, True, False, False])
print(lamarckian_threshold_keys(keys, selected))
# Expected: [0.9 0.9 0.4 0.4] — genes 0 and 2 reflected across 0.5, others kept.

Lamarckian write-back accelerates convergence but reduces diversity; combine it with a higher mutant fraction. The same idea applies after embedding local search in the decoder (a memetic BRKGA): decode, improve with 2-opt/swap descent, write the improved solution's keys back. For online parameter control, the adaptive BRKGA of Chaves, Gonçalves & Lorena (2018, "Adaptive biased random-key genetic algorithm with local search for the capacitated centered clustering problem") shrinks $p_m$ and grows $\rho_e$ over the run, trading early exploration for late intensification — a sensible default schedule when tuning budget is scarce.

Practical Challenges

Decoding dominates runtime and the run crawls. Profile first: in a correct BRKGA, >90% of time is decode. Then (a) vectorize the evaluation across the population as in makespan_batch; (b) exploit the framework's elite-fitness caching — never re-decode elites; (c) parallelize per-chromosome decoders (sequential repairs) across processes; (d) JIT-compile sequential repair loops with numba. Only after these, reduce pop_size.

The population converges to copies of one solution in a few dozen generations. Diversity collapse. Check $p_e \le 0.25$ and $\rho_e \le 0.8$; raise $p_m$ toward 0.3; add duplicate elimination — but compare decoded solutions or fitness fingerprints, not key vectors, because the decoder is many-to-one and distinct keys can hide identical solutions. If collapse recurs, add shaking/restarts.

Best fitness history is not monotone non-increasing. This is a bug, not noise: BRKGA copies elites unchanged, so the generation-best can never get worse. Usual causes: elite fitness not carried over with the elite keys after sorting, a stochastic decoder (decoder must be a deterministic function of the keys — fix any internal randomness with a hash of the chromosome if randomness is truly needed), or fitness recomputed with a different formula than reported.

Warm-start injection makes results worse than cold start. Too many copies of one solution seize the elite set in generation 0 and the run degenerates to mutation-only search around it. Inject 1 exact copy plus 2–3 jittered copies (jittered_copies above), never more than ~5% of the population.

The threshold decoder spends all its time in repair on early generations. With uniform keys, expected initial selection density is 0.5 — far above typical SCP solution density, so pruning does heavy work; in packing-type problems the same bias causes mass infeasibility. Fix the prior: draw initial keys (and mutants) from a Beta distribution or shift the threshold so expected density matches a greedy solution's density.

Sort decoder behaves erratically after key write-back or injection. Tied keys decode by argsort's stable index order, which silently biases sequences. Continuous random keys never tie, but injected or reflected keys can. Add tiny uniform jitter (e.g., $10^{-9}$ scale) when constructing keys from solutions.

BRKGA stalls a few percent above best-known values on routing-flavored problems. Random keys preserve relative order through crossover, not adjacency, so edge-based structure is invisible to the operators. Either embed local search in the decoder (memetic BRKGA, Lamarckian write-back) or switch the representation — this is a model-choice problem, not a tuning problem.

Too many parameters to tune and no budget. Defaults $p = \max(100, n)$, $p_e = 0.2$, $p_m = 0.15$, $\rho_e = 0.7$ are robust starting points (Gonçalves & Resende, 2011). If tuning is worthwhile, tune $\rho_e$ and $p_m$ first; population size mostly trades quality for time. The decoder design moves results far more than any parameter — invest there.

Tools & Libraries

Output Format

A complete BRKGA deliverable contains:

1. Configuration summary — every value needed to reproduce the run:

2. Solution-quality report — final objective, gap to baseline and (if known) to best-known/optimal value, independent validation result. Always state the baseline: "BRKGA 361 vs. greedy 392 (−7.9%); feasibility confirmed by validate_cover."

3. Convergence summary — generation of last improvement, best-per-generation history (plot or table at checkpoints 0/25/50/100%), and whether the stall criterion fired. A run that improves in the final 10% of generations was stopped too early.

4. Runtime breakdown — total wall time, share spent in decode (expected: dominant), decodes per second. This tells the reader whether the next improvement should target the decoder or the search.

5. Replication statistics (when comparing methods) — per-seed final values, mean/median/std and a distribution plot over ≥10 seeds; never report a single lucky run when a comparison claim is made.

6. File artifacts — brkga_core.py (framework), one <problem>.py (instance + evaluator + validator), one run_<problem>.py (decoder + driver), and a results table (CSV) of per-run outcomes.

Questions to Ask

• What does one solution look like structurally — a permutation, a subset, an assignment, or several of these at once?
• Is there an existing constructive heuristic we can drive with priorities (that becomes the decoder)?
• How expensive is evaluating one solution, and what is the total wall-clock budget?
• Which constraints must the decoder absorb or repair, and which can be penalized?
• How large are the instances (chromosome length), and how many instances must be solved?
• What is the quality bar: beat a stated baseline, match published best-known values, or "good solution fast"?
• Are known good solutions available for warm-starting?
• Is a single run enough, or do you need statistics over seeds for a publication-grade comparison?
• Should we use the maintained brkga_mp_ipr package, or is a from-scratch implementation required (custom instrumentation, vectorized decoding)?

Related Skills

• decoder-based-representations — when the decoder itself needs design work: full catalog of sort, threshold, greedy/SGS, allocation, and multi-segment decoders and their feasibility guarantees.
• genetic-algorithms — when BRKGA's fixed operator set is too rigid and you need general GA machinery: selection schemes, permutation crossovers, mutation operators, replacement strategies.
• flow-shop-scheduling — when the application is the flow shop itself: exact MIP models, NEH and other constructive heuristics, lower bounds, and benchmark instances beyond the sort-decoder example here.
• parallel-machine-scheduling — when keys should drive a dispatch/assignment decoder for identical or unrelated parallel machines: list-scheduling decoders, LPT seeds, and makespan/total-completion-time objectives.