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 $ne = \lceil pe \cdot p \rceil$ chromosomes, copied unchanged to the next generation (strong elitism; best fitness is monotone).
• Mutants: $nm = \lceil pm \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 - ne - nm$ 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$:

$$ ci = \begin{cases} ai & \text{with probability } \rhoe \quad (a \in E) \\ bi & \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 $\rhoe = 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 $pe = 0.5` selects, then repair | Set covering, knapsack, facility selection | | Greedy-priority decoder | Any constructive solution | Keys are priorities feeding a constructive heuristic | RCPSP serial SGS, parallel machine dispatch, bin packing | | Allocation decoder | Categorical choice per item | Split $[0,1)$ into intervals, key picks the interval | Machine assignment, mode selection | | Multi-segment chromosome | Several decisions at once | Concatenate key blocks, one decoder stage per block | Assign-then-sequence scheduling, location-then-routing |

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 - ne) \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  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))  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} cj yj \quad \text{s.t.} \quad \sum{j} a{ij}\, yj \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))  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;
    # BR

…

## 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.