Selection And Replacement Strategies

Selection and Replacement Strategies

You are an expert in selection and replacement for evolutionary and population-based metaheuristics. This skill is the reference catalog: for every standard parent-selection scheme — tournament, fitness-proportionate (roulette), linear and exponential rank, stochastic universal sampling, Boltzmann — and every standard survivor strategy — generational, elitist, steady-state, (mu+lambda), (mu,lambda) — it gives when to use it, a numpy implementation, a complexity note, and the algorithms and problem settings it fits. Use the pressure framework below (selection intensity, takeover time) to set pressure deliberately instead of inheriting it by accident from a default operator, and use the measurement harness to verify the setting empirically.

Initial Assessment

Establish these facts before recommending or writing any selection code:

• Maximize or minimize, and can fitness be negative or zero? Fitness-proportionate selection requires strictly positive values and is not invariant to shifting; a penalized objective that dips below zero rules it out unless you accept that the shift itself silently sets the pressure. Rank and tournament are immune.
• Audit every place selection pressure enters. Parent selection, survivor selection, elitism, and duplicate handling all exert pressure, and their effects compound. Tournament k=3 parents feeding a replace-worst steady-state population is a very different machine from the same tournament feeding pure generational replacement.
• Population size N and total evaluation budget. Convert the budget to generations G = budget / lambda, then compare against takeover time. If takeover is far below G, most of the run is post-convergence mutation hill-climbing; if it is far above G, the run never exploits.
• Which algorithm consumes the operator? A canonical GA, a memetic algorithm with strong local search (tolerates and usually wants lower parent pressure), an EDA (wants truncation), an ES with self-adaptive parameters (requires comma replacement), or a multi-objective EA (selection ranks by dominance and crowding, see multi-objective-optimization).
• Generational or steady-state architecture, and how is evaluation parallelized? Generational and (mu,lambda) produce full batches that saturate parallel workers; steady-state evaluates one child at a time and serializes unless you run asynchronous variants.
• Is the objective noisy or stochastic? Noise discounts effective selection intensity and makes naive elitism enshrine lucky evaluations. Plan resampling or comma replacement before tuning anything else.
• What disruption do the variation operators cause? Strong mutation or local search after selection rebuilds diversity each generation and tolerates higher pressure; weak variation under high pressure converges prematurely.
• Any observed pathology? "Converges in 20 generations then stalls" means pressure (or compounded pressure) is too high or diversity machinery is missing; "never improves past random search" usually means pressure is too low or proportionate selection has collapsed.
• Constraint handling. Penalty terms stretch and compress the fitness scale over the run; scale-sensitive schemes (proportionate, raw Boltzmann) inherit that distortion, rank-based schemes do not.
• Reproducibility. Every stochastic scheme takes an explicit np.random.Generator; record seeds per run. A scheme comparison without fixed seeds and repeated runs is noise.

Pressure Anatomy: Intensity, Takeover Time, and the Two Ends of Selection

Two decisions, one budget

Every population algorithm makes two selection decisions per cycle: parent selection (who reproduces, possibly with repeats) and survivor selection / replacement (who stays in the population). Both convert fitness differences into copy-number differences; the product of the two is the algorithm's real pressure. This skill catalogs both ends separately and then measures them together, because tuning one end while ignoring the other is the most common selection mistake.

Selection intensity

The standard scalar for pressure is the selection intensity (Mühlenbein & Schlierkamp-Voosen 1993, "Predictive Models for the Breeder Genetic Algorithm"): the standardized gain in mean fitness produced by selection alone,

$$ I = \frac{\bar f_{\mathrm{sel}} - \bar f}{\sigma_f}, $$

where $\bar f$ and $\sigma_f$ are the population mean and standard deviation and $\bar f_{\mathrm{sel}}$ is the mean of the selected parents. For approximately normal fitness distributions, closed forms exist per scheme (Blickle & Thiele 1996, "A Comparison of Selection Schemes Used in Evolutionary Algorithms") and are tabulated below. Intensity is measurable in any run — log it every generation; it is the cheapest diagnostic in this whole skill.

Takeover time

Takeover time is the number of generations selection alone (no crossover, no mutation) needs before copies of the single best individual fill the population (Goldberg & Deb 1991, "A Comparative Analysis of Selection Schemes Used in Genetic Algorithms"). Approximate forms for population size N:

$$ \tau_{\text{tournament-}k} \approx \frac{\ln N + \ln \ln N}{\ln k}, \qquad \tau_{\text{rank},, s=2} \approx \log_2 N + \log_2 \ln N, \qquad \tau_{\text{proportionate}} \approx N \ln N ;; (f(x) = x), $$

$$ \tau_{(\mu,\lambda)} \approx \frac{\ln \lambda}{\ln(\lambda/\mu)} \qquad \text{(truncation)}. $$

The truncation form is due to Bäck (1996), Evolutionary Algorithms in Theory and Practice. For N = 100: binary tournament takes about 9 generations, tournament k=7 about 3, proportionate selection on linear fitness about 460. The qualitative lesson is robust even where the formulas are loose: ordinal schemes take over in O(log N) generations, proportionate selection in O(N log N), and the gap is the entire design space. Budget rule of thumb: aim for takeover time between G/10 and G/2 of your generation budget G, then verify with the harness below and adjust.

Master catalog: parent selection

Linear rank with s=2 and binary tournament have the same expected behavior (both I = 0.564); tournament is cheaper and needs no global fitness statistics, which is why it became the modern default.

Master catalog: replacement

Scheme-by-context fit

Parent Selection Catalog

All implementations take an explicit np.random.Generator, return integer index arrays into the population (so they compose with any genotype array), and are vectorized — no Python loop touches individuals.

Fitness-proportionate selection (roulette wheel)

When to use. Almost never as a tuning choice in modern practice — its pressure equals the population's coefficient of variation, which you do not control and which decays to zero as the run converges (the "stall" pathology). It survives as the historical baseline (Holland 1975, Adaptation in Natural and Artificial Systems), as the sampling core that rank and Boltzmann reuse, and in settings where fitness genuinely is a nonnegative rate you want copies proportional to. Complexity: O(N) cumulative sum, O(log N) per draw via searchsorted. Fits: teaching, probability-matching components inside EDAs, never penalized objectives.

import numpy as np


def roulette_select(
    fitness: np.ndarray, n_select: int, rng: np.random.Generator
) -> np.ndarray:
    """Fitness-proportionate selection; returns indices of selected parents.

    Requires strictly positive fitness. Each draw is independent, so the
    realized copy counts have high variance; prefer SUS for the same
    expected copies with minimal variance.
    """
    p = fitness / fitness.sum()
    cum = np.cumsum(p)
    cum[-1] = 1.0                                  # guard against rounding
    return np.searchsorted(cum, rng.random(n_select), side="right")


# Tiny instance: fitness 1..5, copies should be proportional to fitness.
rng = np.random.default_rng(0)
fit = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
counts = np.bincount(roulette_select(fit, 15_000, rng), minlength=5)
print(np.round(counts / 15_000, 3))
# Expected: approx [0.067, 0.133, 0.200, 0.267, 0.333] = fitness / 15

Stochastic universal sampling (SUS)

When to use. Whenever you sample parents from an explicit probability vector — proportionate, rank, or Boltzmann probabilities alike. Baker (1987), "Reducing Bias and Inefficiency in the Selection Algorithm": one spin, n_select equally spaced pointers, so every individual receives either the floor or the ceiling of its expected copy count (zero bias, minimal spread). This removes the sampling-variance failure mode of roulette at identical expected pressure. Complexity: O(N) plus one searchsorted over all pointers; effectively free. Fits: drop-in replacement for any wheel.

import numpy as np


def sus_select(
    fitness: np.ndarray, n_select: int, rng: np.random.Generator
) -> np.ndarray:
    """Stochastic universal sampling (Baker 1987) over positive fitness.

    Shuffles the result so downstream pairing of consecutive parents
    is random rather than fitness-sorted.
    """
    p = fitness / fitness.sum()
    cum = np.cumsum(p)
    cum[-1] = 1.0
    pointers = (rng.random() + np.arange(n_select)) / n_select
    return rng.permutation(np.searchsorted(cum, pointers, side="right"))


rng = np.random.default_rng(1)
fit = np.array([1.0, 2.0, 3.0, 4.0])
counts = np.bincount(sus_select(fit, 10, rng), minlength=4)
print(counts)
# Expected: [1 2 3 4] — expected copies n_select * p_i = [1, 2, 3, 4]
# are integers here, so SUS hits them exactly, every time, any seed

Rank selection (linear and exponential)

When to use. When you want proportionate-style smooth probabilities but with pressure that is constant over the run and immune to the fitness scale: rank selection depends only on the ordering, so it is invariant under any monotone transform of the objective (penalties, log-scaling, sign flips handled by sort direction). Linear ranking (Baker 1985, "Adaptive Selection Methods for Genetic Algorithms") gives the best individual s expected copies, s in (1, 2]; the worst gets 2 - s. Exponential ranking breaks the s <= 2 ceiling when you need more pressure but still want smooth probabilities. Complexity: O(N log N) for the sort, then SUS. Fits: GENITOR-style steady-state GAs (Whitley 1989, "The GENITOR Algorithm and Selection Pressure"), any GA whose objective includes penalty terms.

import numpy as np


def linear_rank_probs(fitness: np.ndarray, s: float = 1.7) -> np.ndarray:
    """Baker (1985) linear ranking probabilities; pressure knob 1 < s <= 2.

    Best individual expects s copies per N draws, worst expects 2 - s.
    Depends only on ranks, never on fitness magnitudes.
    """
    n = fitness.size
    ranks = np.empty(n)
    ranks[np.argsort(fitness)] = np.arange(n)      # 0 = worst, n-1 = best
    return (2.0 - s) / n + 2.0 * ranks * (s - 1.0) / (n * (n - 1.0))


def exponential_rank_probs(fitness: np.ndarray, c: float = 0.98) -> np.ndarray:
    """Exponential ranking: weight c^(n-1-rank), c in (0, 1).

    Smaller c means more pressure; c -> 1 approaches uniform selection.
    """
    n = fitness.size
    ranks = np.empty(n, dtype=np.int64)
    ranks[np.argsort(fitness)] = np.arange(n)
    w = c ** (n - 1.0 - ranks)
    return w / w.sum()


def rank_select(
    fitness: np.ndarray, n_select: int, rng: np.random.Generator,
    s: float = 1.7,
) -> np.ndarray:
    """Linear-rank selection drawn with SUS for minimal sampling variance."""
    p = linear_rank_probs(fitness, s)
    cum = np.cumsum(p)
    cum[-1] = 1.0
    pointers = (rng.random() + np.arange(n_select)) / n_select
    return rng.permutation(np.searchsorted(cum, pointers, side="right"))


fit = np.array([10.0, 1e6, 10.5, 11.0])            # one huge outlier
print(np.round(linear_rank_probs(fit, s=2.0), 3))
# Expected: [0. 0.5 0.167 0.333] — the outlier gets rank share 0.5, not
# probability ~1.0 as roulette would give it; ranking absorbs magnitude

Tournament selection

When to use. The modern default. One integer knob k sets pressure (I from 0.56 at k=2 to 1.35 at k=7); no sort, no global fitness statistics, no positivity requirement; it needs only pairwise comparisons, so it works even when "fitness" is the outcome of a head-to-head simulation. The probability that the rank-r individual (r = 0 worst) wins a size-k tournament with replacement is $((r+1)^k - r^k)/N^k$. The probabilistic variant (winner with probability p, else the runner-up, geometrically) interpolates pressure below k=2 when even binary tournament is too greedy. Sokolov & Whitley (2005), "Unbiased Tournament Selection," removes the with-replacement sampling noise by permuting columns instead. Complexity: O(Nk) per generation, embarrassingly parallel. Fits: everything; the first thing to try.

import numpy as np


def tournament_select(
    fitness: np.ndarray, n_select: int, rng: np.random.Generator,
    k: int = 3, p_win: float = 1.0,
) -> np.ndarray:
    """Tournament selection with replacement, fully vectorized.

    Draws an (n_select, k) candidate matrix and takes each row's best.
    With p_win < 1 the rank-j candidate of the tournament wins with
    probability p_win * (1 - p_win)^j (probabilistic tournament).
    """
    cand = rng.integers(0, fitness.size, size=(n_select, k))
    if p_win >= 1.0:
        return cand[np.arange(n_select), np.argmax(fitness[cand], axis=1)]
    order = np.argsort(-fitness[cand], axis=1, kind="stable")
    j = np.minimum(rng.geometric(p_win, size=n_select) - 1, k - 1)
    return cand[np.arange(n_select), order[np.arange(n_select), j]]


rng = np.random.default_rng(3)
fit = np.arange(10, dtype=float)                   # individual i has fitness i
sel = tournament_select(fit, 100_000, rng, k=2)
print(np.round(np.bincount(sel, minlength=10) / 100_000, 3))
# Expected: P(rank r) = ((r+1)^2 - r^2)/100 = (2r+1)/100, i.e. approximately
# [0.01, 0.03, 0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19]

Boltzmann selection

When to use. When you want pressure as an explicit, continuously scheduled quantity: $p_i \propto \exp(f_i / T)$, high temperature T is near-uniform exploration, low T is near-greedy exploitation, and annealing T over the run sweeps the whole range — the population analogue of simulated annealing's acceptance schedule (Goldberg 1990, "A Note on Boltzmann Tournament Selection for Genetic Algorithms and Population-Oriented Simulated Annealing"; Mahfoud & Goldberg 1995, parallel recombinative simulated annealing). Raw Boltzmann is scale-sensitive — dividing fitness by its standard deviation before the exponent makes the temperature scale-free and is strongly recommended. For approximately normal fitness, exponential tilting shifts the mean by $\sigma_f^2/T$, so intensity is about $\sigma_f / T$ (or $1/t$ with the sigma-scaled temperature t). Complexity: O(N). Fits: runs with a known evaluation budget where you want exploration early and exploitation late without changing operators.

import numpy as np


def boltzmann_select(
    fitness: np.ndarray, n_select: int, rng: np.random.Generator,
    temperature: float,
) -> np.ndarray:
    """Boltzmann selection, p_i proportional to exp(f_i / T), computed stably.

    Subtracting the max before exponentiating prevents overflow; sampling
    uses SUS pointers for minimal variance.
    """
    z = (fitness - fitness.max()) / max(temperature, 1e-12)
    p = np.exp(z)
    p /= p.sum()
    cum = np.cumsum(p)
    cum[-1] = 1.0
    pointers = (rng.random() + np.arange(n_select)) / n_select
    return rng.permutation(np.searchsorted(cum, pointers, side="right"))


def geometric_temperature(t0: float, t_end: float, n_steps: int) -> np.ndarray:
    """Geometric annealing schedule from t0 down to t_end over n_steps."""
    alpha = (t_end / t0) ** (1.0 / max(n_steps - 1, 1))
    return t0 * alpha ** np.arange(n_steps)


rng = np.random.default_rng(4)
fit = np.array([1.0, 2.0, 3.0])
for temp in (10.0, 1.0, 0.1):
    counts = np.bincount(boltzmann_select(fit, 9_000, rng, temp), minlength=3)
    print(temp, np.round(counts / 9_000, 2))
# Expected: T=10 near [0.30 0.33 0.37] (almost uniform); T=1 about
# [0.09 0.24 0.67]; T=0.1 about [0.00 0.00 1.00] (effectively greedy)

Measuring pressure directly: intensity and takeover harness

Closed forms assume normal fitness and large N; your penalized, multi-modal, finite-N reality differs. This harness measures both quantities for any scheme: realized intensity from one generation's selection, and takeover time by running selection alone (no variation) until the best type holds 90% of the population. The 90% criterion (instead of 100%) keeps near-neutral end-game drift from dominating the statistic; the best starts with 5% of the slots because a single copy goes extinct with noticeable probability under high-variance schemes, which would bias the median upward.

import numpy as np
from collections.abc import Callable

SelectFn = Callable[[np.ndarray, int, np.random.Generator], np.ndarray]


def realized_intensity(fitness: np.ndarray, selected: np.ndarray) -> float:
    """Realized selection intensity I = (mean(selected) - mean) / std."""
    sd = float(fitness.std())
    if sd == 0.0:
        return 0.0
    return float((fitness[selected].mean() - fitness.mean()) / sd)


def takeover_time(
    select: SelectFn, n_pop: int, rng: np.random.Generator,
    share_target: float = 0.9, max_gen: int = 5_000,
) -> int:
    """Generations of selection alone until the best type holds share_target.

    Selection-only dynamics in the sense of Goldberg & Deb (1991): the
    population is just a fitness multiset, each generation replaced by the
    selected multiset. Returns max_gen if the target share is never reached
    (e.g. the best type drifted to extinction under a high-variance scheme).
    """
    fitness = np.linspace(0.5, 0.99, n_pop)
    fitness[: max(1, n_pop // 20)] = 1.0           # best starts at 5% share
    for gen in range(1, max_gen + 1):
        fitness = fitness[select(fitness, n_pop, rng)]
        if (fitness == 1.0).mean() >= share_target:
            return gen
    return max_gen


def median_takeover(select: SelectFn, n_pop: int, n_runs: int, seed: int) -> float:
    """Median takeover time over independent replicates."""
    times = [takeover_time(select, n_pop, np.random.default_rng(seed + r))
             for r in range(n_runs)]
    return float(np.median(times))


def tournament_k(k: int) -> SelectFn:
    """Size-k tournament closure for the harness."""
    def sel(f: np.ndarray, n: int, rng: np.random.Generator) -> np.ndarray:
        cand = rng.integers(0, f.size, size=(n, k))
        return cand[np.arange(n), np.argmax(f[cand], axis=1)]
    return sel


def proportionate(f: np.ndarray, n: int, rng: np.random.Generator) -> np.ndarray:
    """Roulette wheel for the harness (fitness here is already positive)."""
    cum = np.cumsum(f / f.sum())
    cum[-1] = 1.0
    return np.searchsorted(cum, rng.random(n), side="right")


for name, fn in [("tournament k=2", tournament_k(2)),
                 ("tournament k=7", tournament_k(7)),
                 ("proportionate", proportionate)]:
    print(name, median_takeover(fn, n_pop=100, n_runs=9, seed=7))
# Expected: tournament k=2 about 6-9 generations and k=7 about 2-4
# (Goldberg-Deb predict 8.9 and 3.2 for full takeover from a single copy at
# N=100); proportionate is one to two orders of magnitude slower — tens to a
# few hundred generations — because its pressure equals the population's
# relative fitness spread, which shrinks exactly as takeover proceeds.

Replacement and Survivor Selection Catalog

Replacement decides the generation gap (De Jong 1975 dissertation: how much of the population turns over per cycle) and whether the incumbent can be lost. Generational replacement swaps the whole population and relies entirely on parent selection for pressure; adding elitism (copy the e best parents over the worst offspring) is the cheap fix that makes best-so-far monotone inside the population. Steady-state replacement (Whitley 1989, GENITOR; Syswerda 1991, "A Study of Reproduction in Generational and Steady-State Genetic Algorithms") inserts one child at a time: replace-worst adds strong survivor pressure on top of parent pressure, while replace-random or replace-oldest behaves close to generational. (mu+lambda) pools parents and offspring and keeps the best mu — pure exploitation, never loses the incumbent. (mu,lambda) keeps the best mu of the lambda offspring only — it deliberately forgets, which is exactly what self-adaptive strategy parameters and noisy or time-varying objectives need (Beyer & Schwefel 2002, "Evolution Strategies — A Comprehensive Introduction"); keep an external best-so-far archive for reporting. Truncation pressure is set by lambda/mu (intensity 0.80 at tau=1/2, 1.27 at 1/4, 1.75 at 1/10).

Complexity: all of the following are O(N) or O(N log N) per generation with argpartition; replacement is never the bottleneck. Fits: see the replacement table in the anatomy section.

import numpy as np


def generational_replace(
    pop: np.ndarray, fit: np.ndarray,
    off: np.ndarray, off_fit: np.ndarray, n_elite: int = 1,
) -> tuple[np.ndarray, np.ndarray]:
    """Offspring replace the population; the n_elite best parents survive
    by overwriting the worst offspring (elitism, De Jong 1975)."""
    new_pop, new_fit = off.copy(), off_fit.copy()
    if n_elite > 0:
        elite = np.argpartition(-fit, n_elite - 1)[:n_elite]
        victims = np.argpartition(new_fit, n_elite - 1)[:n_elite]
        new_pop[victims], new_fit[victims] = pop[elite], fit[elite]
    return new_pop, new_fit


def mu_plus_lambda_replace(
    pop: np.ndarray, fit: np.ndarray,
    off: np.ndarray, off_fit: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """(mu+lambda): keep the best mu = len(pop) of parents and offspring pooled."""
    mu = pop.shape[0]
    all_pop = np.vstack((pop, off))
    all_fit = np.concatenate((fit, off_fit))
    keep = np.argpartition(-all_fit, mu - 1)[:mu]
    return all_pop[keep], all_fit[keep]


def mu_comma_lambda_replace(
    off: np.ndarray, off_fit: np.ndarray, mu: int,
) -> tuple[np.ndarray, np.ndarray]:
    """(mu,lambda): keep the best mu offspring; all parents die.

    Requires lambda >= mu (ES practice: lambda/mu in 4-7). Can lose the
    best-so-far by design — archive the incumbent externally.
    """
    keep = np.argpartition(-off_fit, mu - 1)[:mu]
    return off[keep], off_fit[keep]


def steady_state_insert(
    pop: np.ndarray, fit: np.ndarray, child: np.ndarray, child_fit: float,
) -> bool:
    """Replace-worst steady-state insertion with duplicate rejection.

    Rejects genotypes already present — without this check replace-worst
    floods the population with copies. Modifies pop/fit in place; returns
    True if the child was inserted.
    """
    if bool((pop == child).all(axis=1).any()):
        return False
    worst = int(np.argmin(fit))
    if child_fit <= fit[worst]:
        return False
    pop[worst], fit[worst] = child, child_fit
    return True


pop = np.array([[1, 1, 0], [0, 1, 1], [1, 0, 0], [0, 0, 0]], dtype=np.int8)
fit = np.array([8.0, 7.0, 3.0, 1.0])
off = np.array([[1, 1, 1], [0, 1, 0], [1, 0, 1], [0, 0, 1]], dtype=np.int8)
off_fit = np.array([9.0, 2.0, 6.0, 4.0])

print(sorted(generational_replace(pop, fit, off, off_fit, n_elite=1)[1]))
# Expected: [4.0, 6.0, 8.0, 9.0] — elite parent 8.0 replaced worst child 2.0
print(sorted(mu_plus_lambda_replace(pop, fit, off, off_fit)[1]))
# Expected: [6.0, 7.0, 8.0, 9.0] — best four of the pooled eight
print(sorted(mu_comma_lambda_replace(off, off_fit, mu=2)[1]))
# Expected: [6.0, 9.0] — best two offspring; parent 8.0 is forgotten
print(steady_state_insert(pop, fit, np.array([1, 1, 0], dtype=np.int8), 9.0))
# Expected: False — genotype already in the population (duplicate rejection)
print(steady_state_insert(pop, fit, np.array([1, 1, 1], dtype=np.int8), 9.0))
# Expected: True — replaces the worst (fitness 1.0) in place

Worked Example: Selection-Pressure Experiment on One GA

The catalog's promise is that intensity predicts behavior. This experiment checks it on one fixed GA — generational with 1 elite, uniform crossover, 1/n bit-flip mutation — solving a penalized 0-1 knapsack, swapping only the parent-selection scheme. Two outputs per scheme: the mean realized intensity over the run (the measured knob) and the best objective over seeds (the consequence). Proportionate selection uses SUS sampling and a min-shift, which is the most favorable implementation it can get; it still loses because the shift plus the shrinking fitness spread destroy its pressure late in the run. Boltzmann uses a sigma-scaled geometric temperature schedule, so its intensity ramps from about 0.5 to effectively greedy.

import numpy as np
import pandas as pd
from collections.abc import Callable

SelectFn = Callable[[np.ndarray, int, np.random.Generator, int], np.ndarray]


def make_knapsack(n_items: int, seed: int) -> tuple[np.ndarray, np.ndarray, float]:
    """Weakly correlated 0-1 knapsack: values tied to weights, capacity 50%."""
    rng = np.random.default_rng(seed)
    weights = rng.integers(5, 30, n_items).astype(float)
    values = weights + rng.integers(1, 15, n_items)
    return values, weights, 0.5 * float(weights.sum())


def fitness_pop(
    pop: np.ndarray, values: np.ndarray, weights: np.ndarray, capacity: float
) -> np.ndarray:
    """Penalized fitness of an (N, n) 0/1 population.

    The penalty rate exceeds every item's value/weight ratio, so no
    infeasible solution can outscore the feasible optimum.
    """
    rho = float((values / weights).max()) + 1.0
    return pop @ values - rho * np.maximum(pop @ weights - capacity, 0.0)


def sus_from_probs(p: np.ndarray, n: int, rng: np.random.Generator) -> np.ndarray:
    """SUS draw of n indices from a probability vector."""
    cum = np.cumsum(p)
    cum[-1] = 1.0
    idx = np.searchsorted(cum, (rng.random() + np.arange(n)) / n, side="right")
    return rng.permutation(idx)


def make_schemes(n_gens: int) -> dict[str, SelectFn]:
    """Five schemes under test, all with signature (fitness, n, rng, gen)."""

    def prop(f: np.ndarray, n: int, rng: np.random.Generator,
             gen: int) -> np.ndarray:
        shifted = f - f.min() + 1e-9               # proportionate needs > 0
        return sus_from_probs(shifted / shifted.sum(), n, rng)

    def lin_rank(f: np.ndarray, n: int, rng: np.random.Generator,
                 gen: int, s: float = 1.7) -> np.ndarray:
        m = f.size
        ranks = np.empty(m)
        ranks[np.argsort(f)] = np.arange(m)
        p = (2.0 - s) / m + 2.0 * ranks * (s - 1.0) / (m * (m - 1.0))
        return sus_from_probs(p, n, rng)

    def tournament(k: int) -> SelectFn:
        def sel(f: np.ndarray, n: int, rng: np.random.Generator,
                gen: int) -> np.ndarray:
            cand = rng.integers(0, f.size, size=(n, k))
            return cand[np.arange(n), np.argmax(f[cand], axis=1)]
        return sel

    def boltzmann(f: np.ndarray, n: int, rng: np.random.Generator,
                  gen: int) -> np.ndarray:
        t = 2.0 * (0.02 / 2.0) ** (gen / max(n_gens - 1, 1))   # 2.0 -> 0.02
        z = (f - f.max()) / (t * max(float(f.std()), 1e-9))    # sigma-scaled
        p = np.exp(z)
        return sus_from_probs(p / p.sum(), n, rng)

    return {"proportionate (SUS)": prop, "linear rank s=1.7": lin_rank,
            "tournament k=2": tournament(2), "tournament k=7": tournament(7),
            "Boltzmann t 2->0.02": boltzmann}


def run_ga(
    select: SelectFn, values: np.ndarray, weights: np.ndarray, capacity: float,
    n_pop: int = 100, n_gens: int = 150, seed: int = 0,
) -> tuple[float, float]:
    """Generational GA (uniform crossover, 1/n bit-flip, 1 elite).

    Returns (best penalized objective, mean realized selection intensity).
    """
    rng = np.random.default_rng(seed)
    n = values.size
    pop = (rng.random((n_pop, n)) < 0.5).astype(np.int8)
    fit = fitness_pop(pop, values, weights, capacity)
    best_i = int(np.argmax(fit))
    best_x, best_f = pop[best_i].copy(), float(fit[best_i])
    intensity = np.empty(n_gens)
    for gen in range(n_gens):
        parents = select(fit, n_pop, rng, gen)
        sd = float(fit.std())
        intensity[gen] = 0.0 if sd == 0 else (fit[parents].mean() - fit.mean()) / sd
        pa, pb = pop[parents[0::2]], pop[parents[1::2]]
        mask = rng.random(pa.shape) < 0.5
        off = np.vstack((np.where(mask, pa, pb), np.where(mask, pb, pa)))
        off ^= (rng.random(off.shape) < 1.0 / n).astype(np.int8)
        fit = fitness_pop(off, values, weights, capacity)
        victim = int(np.argmin(fit))               # 1-elite replacement
        off[victim], fit[victim] = best_x, best_f
        pop = off
        g = int(np.argmax(fit))
        if fit[g] > best_f:
            best_x, best_f = pop[g].copy(), float(fit[g])
    return best_f, float(intensity.mean())


values, weights, capacity = make_knapsack(60, seed=42)
rows = []
for name, sel in make_schemes(n_gens=150).items():
    out = np.array([run_ga(sel, values, weights, capacity, seed=s)
                    for s in range(5)])
    rows.append({"scheme": name,
                 "mean_best": round(float(out[:, 0].mean()), 1),
                 "std_best": round(float(out[:, 0].std()), 1),
                 "mean_intensity": round(float(out[:, 1].mean()), 2)})
print(pd.DataFrame(rows).to_string(index=False))
# Expected: measured mean_intensity orders the schemes
# proportionate < rank s=1.7 < tournament k=2 < tournament k=7 (about
# 0.33 / 0.39 / 0.55 / 0.96 here), with Boltzmann's per-generation intensity
# sweeping from ~0.5 to greedy as t anneals. mean_best follows the same
# order: proportionate trails with the largest std_best — its pressure
# collapses once the population's relative fitness spread narrows — while
# the high-pressure schemes win on this small, easy instance. On deceptive
# or rugged landscapes the high-pressure end pays for that speed with
# premature convergence; verify with the takeover harness before copying k=7.

Read the table the way you would tune a real run: if mean_intensity is below about 0.3 the scheme is barely selecting (raise s, k, or lower t); if it is above about 1.3 on a problem this small, expect seed-to-seed variance and premature convergence, and check the diversity diagnostics from diversity-and-population-management before trusting the best run.

Worked Example: Replacement Strategies with Convergence Plots

Now hold parent selection fixed (tournament k=3) and swap only the survivor strategy, comparing at equal evaluation budget — the only fair x-axis, since (mu,lambda) here spends 2N evaluations per generation while steady-state spends 1 per step. The plot shows mean best-so-far over seeds with a min-max band per strategy.

import numpy as np
import matplotlib.pyplot as plt


def make_knapsack(n_items: int, seed: int) -> tuple[np.ndarray, np.ndarray, float]:
    """Weakly correlated 0-1 knapsack instance."""
    rng = np.random.default_rng(seed)
    weights = rng.integers(5, 30, n_items).astype(float)
    values = weights + rng.integers(1, 15, n_items)
    return values, weights, 0.5 * float(weights.sum())


def fitness_pop(pop: np.ndarray, values: np.ndarray, weights: np.ndarray,
                capacity: float) -> np.ndarray:
    """Penalized fitness of an (N, n) 0/1 population."""
    rho = float((values / weights).max()) + 1.0
    return pop @ values - rho * np.maximum(pop @ weights - capacity, 0.0)


def make_offspring(pop: np.ndarray, fit: np.ndarray, n_off: int,
                   rng: np.random.Generator) -> np.ndarray:
    """n_off children: tournament k=3 parents, uniform crossover, 1/n bit-flip."""
    cand = rng.integers(0, fit.size, size=(2 * n_off, 3))
    parents = cand[np.arange(2 * n_off), np.argmax(fit[cand], axis=1)]
    pa, pb = pop[parents[:n_off]], pop[parents[n_off:]]
    child = np.where(rng.random(pa.shape) < 0.5, pa, pb)
    return child ^ (rng.random(child.shape) < 1.0 / pop.shape[1]).astype(np.int8)


def run(strategy: str, values: np.ndarray, weights: np.ndarray, capacity: float,
        seed: int, n_pop: int = 100, budget: int = 20_000,
        ) -> tuple[np.ndarray, np.ndarray]:
    """One GA run under a replacement strategy; returns (evals, best_so_far)."""
    rng = np.random.default_rng(seed)
    pop = (rng.random((n_pop, values.size)) < 0.5).astype(np.int8)
    fit = fitness_pop(pop, values, weights, capacity)
    evals, best = n_pop, float(fit.max())
    ev_log, best_log = [float(evals)], [best]
    while evals < budget:
        if strategy == "steady-state":
            child = make_offspring(pop, fit, 1, rng)[0]
            cf = float(fitness_pop(child[None, :], values, weights, capacity)[0])
            evals += 1
            if cf > fit.min() and not bool((pop == child).all(axis=1).any()):
                w = int(np.argmin(fit))
                pop[w], fit[w] = child, cf
        elif strategy == "mu,lambda":
            off = make_offspring(pop, fit, 2 * n_pop, rng)
            of = fitness_pop(off, values, weights, capacity)
            evals += 2 * n_pop
            keep = np.argpartition(-of, n_pop - 1)[:n_pop]
            pop, fit = off[keep], of[keep]
        else:
            off = make_offspring(pop, fit, n_pop, rng)
            of = fitness_pop(off, values, weights, capacity)
            evals += n_pop
            if strategy == "mu+lambda":
                allp = np.vstack((pop, off))
                allf = np.concatenate((fit, of))
                keep = np.argpartition(-allf, n_pop - 1)[:n_pop]
                pop, fit = allp[keep], allf[keep]
            elif strategy == "generational+1elite":
                v, e = int(np.argmin(of)), int(np.argmax(fit))
                off[v], of[v] = pop[e], fit[e]
                pop, fit = off, of
            else:                                   # pure generational
                pop, fit = off, of
        best = max(best, float(fit.max()))
        if evals % n_pop == 0 or evals >= budget:
            ev_log.append(float(evals))
            best_log.append(best)
    return np.array(ev_log), np.array(best_log)


values, weights, capacity = make_knapsack(60, seed=42)
strategies = ["generational", "generational+1elite", "steady-state",
              "mu+lambda", "mu,lambda"]
grid = np.arange(100, 20_001, 100, dtype=float)
fig, ax = plt.subplots(figsize=(7.0, 4.2))
summary: dict[str, float] = {}
for strat in strategies:
    curves = np.vstack([
        np.interp(grid, *run(strat, values, weights, capacity, seed=s))
        for s in range(5)])
    mean = curves.mean(axis=0)
    ax.plot(grid, mean, label=strat)
    ax.fill_between(grid, curves.min(axis=0), curves.max(axis=0), alpha=0.15)
    summary[strat] = round(float(mean[-1]), 1)
ax.set_xlabel("objective evaluations")
ax.set_ylabel("best penalized objective so far")
ax.legend(frameon=False)
fig.tight_layout()
fig.savefig("replacement_comparison.png", dpi=200)
print(summary)
# Expected: steady-state leads from the first checkpoints (it recycles good
# survivors after every single evaluation) and finishes best; mu+lambda is
# next fastest early but plateaus soonest — pure exploitation converges
# first; mu,lambda starts slower and closes the gap late, since its doubled
# lambda per generation buys exploration; pure generational trails the whole
# way because the best is repeatedly lost and rediscovered. The four
# pressure-bearing strategies end within ~0.3% of each other: on a static,
# noise-free instance of this size, replacement choice moves convergence
# speed far more than final quality.

Do not over-read a single instance: the ordering among the pressure-bearing strategies is budget- and problem-dependent, and on noisy objectives the comparison changes character entirely (naive elitism and the plus strategy lock in lucky evaluations; comma flushes them). The reusable artifact is the harness — equal-evaluation accounting, multiple seeds, best-so-far curves with bands — not this instance's specific ranking.

Advanced Techniques

Fitness scaling for proportionate selection

If you must keep a proportionate wheel (legacy code, probability-matching semantics), control its pressure with scaling (Goldberg 1989, Genetic Algorithms in Search, Optimization and Machine Learning, ch. 4). Linear scaling sets $f' = af + b$ with a, b chosen each generation so $\bar f' = \bar f$ and $f'_{\max} = C \bar f$, pinning the best individual's expected copies to C (typical C in 1.2-2.0, cl

…(truncated)