Constraint Handling Techniques

Constraint-Handling Techniques

You are an expert in constraint handling for metaheuristics and evolutionary computation. This skill catalogs the six main technique families — penalty functions (static, dynamic, adaptive), repair operators, feasibility-preserving operators, decoder-based feasibility, stochastic ranking, and Deb's feasibility rules — with numpy implementations, complexity notes, and per-constraint-type selection guidance. Use the framework below to pick a technique per constraint, implement it correctly, and verify the choice empirically with a head-to-head experiment.

Initial Assessment

Establish these facts before recommending any technique:

• Constraint inventory. List every constraint. For each: inequality or equality? Linear or black-box? How many?
• Hard vs soft. Hard constraints define feasibility; soft constraints are preferences. Soft constraints belong in the objective (weighted or lexicographic), never in a feasibility mechanism. Confirm the user agrees on the split.
• Feasible-region density. Sample random solutions: what fraction is feasible? Above ~10%, penalties and feasibility rules work out of the box. Below ~0.1%, you need repair, decoders, or feasibility-preserving operators — random search will never find the feasible region.
• Constraint structure. Is feasibility cheap to check (O(n) capacity sums) or expensive (a simulation)? Cheap checks enable repair and move filtering; expensive checks favor penalties on cached violation values.
• Representation already chosen? If the encoding is still open, the cheapest fix is to encode constraints away (permutation encoding for "visit each once", fixed-cardinality sets for "choose exactly k"). See solution-encodings before adding machinery here.
• Algorithm family. Population methods (GA, DE, EDA) can rank by violation across a population; single-solution methods (SA, tabu, ILS) need per-move decisions — repair, move filtering, or a penalized delta.
• Equality constraints present? Penalties handle equalities poorly (the feasible set has measure zero). Plan for reformulation, decoders, or projection-style repair.
• Where does the optimum live? For most resource-constrained problems the optimum sits on the feasibility boundary (Michalewicz & Schoenauer 1996). Techniques that cannot search near or across the boundary lose quality.
• Evaluation budget and time limit. Repair and decoding add per-individual cost; confirm the budget tolerates it.
• Validation hook. Confirm an independent feasibility checker exists (separate from the fitness code) so the chosen technique can be audited — see solution-validation-testing.

Technique Taxonomy and Selection

The constrained combinatorial problem, in minimization form:

$$ \min{x \in S} f(x) \quad \text{s.t.} \quad gj(x) \le 0 \;(j = 1,\dots,m), \qquad h_k(x) = 0 \;(k = 1,\dots,p), $$

where $S$ is the space reachable by the representation. Define per-constraint violations

$$ vj(x) = \max(0,\, gj(x)), \qquad v{m+k}(x) = \max(0,\, |hk(x)| - \varepsilon), $$

with a tolerance $\varepsilon > 0$ for equalities, and the scalar total violation

$$ \phi(x) = \sumj \left( \frac{vj(x)}{s_j} \right)^{\beta}, $$

where $s_j$ is a per-constraint scale (typical: the constraint's right-hand side or a sampled violation quantile) and $\beta \in \{1, 2\}$. Every technique below is a different way to use $f$ and $\phi$ (or to avoid ever creating $\phi > 0$).

Six families (Coello 2002, "Theoretical and numerical constraint-handling techniques used with evolutionary algorithms"; Mezura-Montes & Coello 2011 survey):

1. Penalty functions — optimize $F(x) = f(x) + \rho \cdot \phi(x)$. Static $\rho$, time-dependent $\rho(t)$, or population-adaptive $\rho$.
2. Feasibility rules (Deb 2000) — never aggregate: compare solutions lexicographically by $(\phi, f)$.
3. Stochastic ranking (Runarsson & Yao 2000) — probabilistic balance between $f$-comparisons and $\phi$-comparisons during sorting.
4. Repair operators — map infeasible solutions to feasible ones after variation.
5. Decoder-based feasibility — the genotype-to-phenotype mapping only produces feasible solutions.
6. Feasibility-preserving operators — variation moves are filtered or designed so feasible parents always yield feasible offspring.

Technique comparison

| Technique | Tuning burden | Offspring always feasible | Searches infeasible region | Extra cost per evaluation | Main failure mode | |---|---|---|---|---|---| | Static penalty | High ($\rho$, $\beta$) | No | Yes | $O(m)$ | $\rho$ too low → infeasible winner; too high → stuck at boundary | | Dynamic penalty | Medium ($C$, $\alpha$) | No (late: nearly) | Yes, early | $O(m)$ | Schedule mismatched to evaluation budget | | Adaptive penalty (APM, Bean–Hadj-Alouane) | Low to none | No | Yes | $O(\text{pop} \cdot m)$ per generation | Coefficient oscillation on noisy populations | | Deb feasibility rules | None | No | Barely | $O(m)$ | Diversity collapse onto first feasible basin | | Stochastic ranking | Low ($Pf \approx 0.45$) | No | Controlled | $O(\text{pop} \cdot \text{sweeps})$ per ranking | Sensitive $Pf$ on very tight problems | | Repair | Design effort | Yes (after repair) | Variation yes, evaluation no | Repair routine cost | Repair bias: offspring cluster on repair targets | | Decoder | Design effort | Yes | No | Decode cost per individual | Locality loss, genotype redundancy | | Feasibility-preserving operators | Design effort | Yes | No | Move feasibility check | Feasible region disconnected under the move set |

Choosing per constraint type

| Constraint type | Example | First choice | Second choice | |---|---|---|---| | Representation-level (each-exactly-once, cardinality $=k$) | TSP tour, choose $k$ medians | Encode it away (permutation / fixed-size set) | Feasibility-preserving operators | | Single resource budget | 0-1 knapsack capacity | Greedy repair (Chu & Beasley 1998) | Calibrated static penalty | | Several resource budgets | Multidimensional knapsack, GAP | Repair + Deb rules for ranking | Adaptive penalty (APM) | | Coupling / equality | Flow balance, demand = supply | Decoder or reformulate variables out | Projection-style repair | | Sparse feasible region ( np.ndarray: """Elementwise inequality violations max(0, g_j(x)) for a batch.

G: (pop, m) constraint values with the convention feasible gj np.ndarray: """Scalar violation per individual: sumj (vj / sj)^beta.""" V = np.maximum(G, 0.0) if scale is not None: V = V / scale return (V**beta).sum(axis=1)

Tiny instance: 3 individuals, 2 constraints (g np.ndarray:

"""Penalized objective F = f + rho sumj max(0, gj)^beta (minimization).""" V = np.maximum(G, 0.0) return f + rho (V**beta).sum(axis=1)

def calibraterho( fsample: np.ndarray, G_sample: np.ndarray, margin: float = 2.0 ) -> float: """Set rho from a random sample so typical violations outweigh the f range.

Rule: rho = margin (objective range) / (mean violation of infeasible samples). """ V = np.maximum(Gsample, 0.0).sum(axis=1) infeasible = V > 0 if not infeasible.any(): return 1.0 frange = float(fsample.max() - fsample.min()) + 1e-12 return margin f_range / float(V[infeasible].mean())

f = np.array([100.0, 90.0, 80.0]) G = np.array([[0.0, 0.0], [2.0, 0.0], [3.0, 1.0]]) print(static_penalty(f, G, rho=10.0))

Expected: [100. 130. 180.] -- the feasible solution now ranks first

**Complexity.** $O(\text{pop} \cdot m)$ per generation. **Fits.** Any algorithm; the only option when the algorithm needs one scalar fitness (e.g., roulette selection, SA acceptance on deltas). Theory note: for linearly constrained problems an exact penalty exists once $\rho$ exceeds the largest dual multiplier, but in metaheuristics $\rho$ is a search-control parameter, not just a correctness parameter — too-large $\rho$ makes all infeasible solutions equally hopeless and freezes the search at the feasibility boundary.

### Dynamic and annealing penalties

**When to use.** You want early exploration through infeasible space and late feasibility pressure, and you know the total generation budget (the schedule must be tied to it).

```python
import numpy as np

def dynamic_penalty(
    f: np.ndarray,
    G: np.ndarray,
    t: int,
    C: float = 0.5,
    alpha: float = 2.0,
    beta: float = 2.0,
) -> np.ndarray:
    """Joines & Houck (1994): F = f + (C*t)^alpha * sum_j v_j^beta, t = generation >= 1."""
    V = np.maximum(G, 0.0)
    return f + (C * t) ** alpha * (V**beta).sum(axis=1)

def annealing_penalty(
    f: np.ndarray, G: np.ndarray, tau: float, beta: float = 2.0
) -> np.ndarray:
    """Michalewicz & Attia (1994): F = f + sum_j v_j^beta / (2*tau), tau cooled toward 0."""
    V = np.maximum(G, 0.0)
    return f + (V**beta).sum(axis=1) / (2.0 * tau)

f = np.array([10.0, 9.0])  # second solution is better on f ...
G = np.array([[0.0], [1.0]])  # ... but violates the constraint by 1
for t in (1, 10, 50):
    print(t, dynamic_penalty(f, G, t))
# Expected: at t=1 the infeasible point still wins (9.25  10.0) and the feasible point wins

Complexity. Same as static per generation. Fits. Population methods with a fixed generation budget; SA (annealing penalty pairs naturally with the SA temperature). Pitfall: with $(Ct)^{\alpha}$ and $\alpha = 2$, pressure grows quadratically — if the budget doubles, the schedule must be re-tuned or the run spends most generations in the effectively-static high-penalty regime.

Adaptive penalties: Bean–Hadj-Alouane and APM

When to use. You cannot calibrate $\rho$ offline, instance scales vary, or feasibility difficulty changes during the run. Adaptive schemes read the population and set coefficients automatically.

import numpy as np

class PenaltyController:
    """Bean & Hadj-Alouane (1992): multiply/divide rho from recent best-individual feasibility.

    If the generation-best was feasible for `window` straight generations, rho /= factor
    (relax, search the boundary); if infeasible for `window` straight, rho *= factor.
    """

    def __init__(self, rho: float = 1.0, factor: float = 2.0, window: int = 5):
        self.rho = rho
        self.factor = factor
        self.window = window
        self.history: list[bool] = []

    def update(self, best_is_feasible: bool) -> float:
        """Record this generation's best-individual feasibility; return updated rho."""
        self.history.append(best_is_feasible)
        recent = self.history[-self.window :]
        if len(recent) == self.window:
            if all(recent):
                self.rho /= self.factor
            elif not any(recent):
                self.rho *= self.factor
        return self.rho

def apm_fitness(f: np.ndarray, G: np.ndarray) -> np.ndarray:
    """Adaptive Penalty Method, Barbosa & Lemonge (2003): parameter-free coefficients.

    k_j = |mean f| * mean(v_j) / sum_l mean(v_l)^2, recomputed every generation.
    Infeasible solutions with f below the population mean are first lifted to the mean.
    """
    V = np.maximum(G, 0.0)
    f_mean = float(f.mean())
    v_mean = V.mean(axis=0)
    denom = float((v_mean**2).sum())
    if denom == 0.0:  # fully feasible population: no penalty needed
        return f.copy()
    k = abs(f_mean) * v_mean / denom
    lifted = np.where(f > f_mean, f, f_mean)
    feasible = V.sum(axis=1) == 0
    return np.where(feasible, f, lifted + V @ k)

f = np.array([5.0, 3.0, 7.0, 4.0])
G = np.array([[0.0], [2.0], [0.0], [1.0]])
print(np.round(apm_fitness(f, G), 2))
# Expected: [ 5.   17.42  7.   11.08] -- feasible kept, infeasible pushed past mean f

Complexity. APM adds $O(\text{pop} \cdot m)$ statistics per generation. Fits. GAs and DE on problems with several heterogeneous constraints; APM also self-scales across constraints, which removes the manual $s_j$ choice. Pitfall: both schemes assume the population statistics are meaningful — with tiny populations (under ~20) the coefficients oscillate; smooth them with an exponential moving average.

Ranking, Repair, and Construction-Based Techniques

Deb's feasibility rules

When to use. Default parameter-free choice for population methods. Three rules (Deb 2000, "An efficient constraint handling method for genetic algorithms"): feasible beats infeasible; two feasible compare on $f$; two infeasible compare on total violation $\phi$. Equivalent to lexicographic order on $(\phi, f)$ — no penalty coefficient exists at all.

import numpy as np

def deb_rank(f: np.ndarray, v: np.ndarray) -> np.ndarray:
    """Indices best-first under Deb (2000): sort by (total violation, objective)."""
    return np.lexsort((f, v))

def deb_tournament(
    f: np.ndarray, v: np.ndarray, n_offspring: int, rng: np.random.Generator
) -> np.ndarray:
    """Vectorized binary tournament under the feasibility rules; returns winner indices."""
    n = len(f)
    a = rng.integers(n, size=n_offspring)
    b = rng.integers(n, size=n_offspring)
    a_wins = (v[a]  np.ndarray:
    """Runarsson & Yao (2000): index array best-first via stochastic bubble sort.

    The sort is inherently sequential; sweeps default to the population size.
    """
    rng = np.random.default_rng() if rng is None else rng
    n = len(f)
    sweeps = n if sweeps is None else sweeps
    idx = np.arange(n)
    for _ in range(sweeps):
        u = rng.random(n - 1)
        swapped = False
        for i in range(n - 1):
            a, b = idx[i], idx[i + 1]
            both_feasible = v[a] == 0.0 and v[b] == 0.0
            compare_on_f = both_feasible or u[i]  f[b]) if compare_on_f else (v[a] > v[b])
            if out_of_order:
                idx[i], idx[i + 1] = b, a
                swapped = True
        if not swapped:
            break
    return idx

f = np.array([1.0, 2.0, 3.0, 4.0])
v = np.array([5.0, 0.0, 0.0, 1.0])
rng = np.random.default_rng(0)
print(stochastic_ranking(f, v, p_f=0.0, rng=rng))  # pure feasibility-first
print(stochastic_ranking(f, v, p_f=1.0, rng=rng))  # pure objective order
# Expected: [1 2 3 0] for p_f=0.0 and [0 1 2 3] for p_f=1.0; p_f=0.45 interpolates

Complexity. $O(\text{pop}^2)$ worst case per ranking (bubble sort), fine for populations up to a few thousand. Fits. $(\mu, \lambda)$ evolution strategies (its original home), GAs with rank-based selection. Tuning: $Pf \in [0.4, 0.475]$; $Pf \ge 0.5$ loses the feasibility guarantee in expectation.

Repair operators

When to use. A cheap greedy map from infeasible to feasible exists. For budget constraints this is the empirical winner — on multidimensional knapsack, repair-based GAs dominate penalty GAs (Chu & Beasley 1998, "A genetic algorithm for the multidimensional knapsack problem").

import numpy as np

def repair_knapsack(
    X: np.ndarray, w: np.ndarray, p: np.ndarray, cap: float
) -> np.ndarray:
    """Greedy DROP/ADD repair for 0-1 knapsack populations (Chu & Beasley 1998).

    DROP: remove worst profit/weight items until the load fits.
    ADD:  insert best-ratio missing items that still fit (also improves feasible rows).
    Returns a repaired copy; X is (pop, n) with 0/1 entries.
    """
    X = X.copy()
    ratio = p / w
    drop_order = np.argsort(ratio)  # worst ratio first
    add_order = drop_order[::-1]  # best ratio first
    load = X @ w
    for i in np.flatnonzero(load > cap):
        for j in drop_order:
            if X[i, j]:
                X[i, j] = 0
                load[i] -= w[j]
                if load[i]  np.ndarray:
    """Random-key decoder for 0-1 knapsack: every genotype maps to a feasible solution.

    Items are considered in decreasing key order and inserted while they fit.
    K: (pop, n) real-valued keys in [0, 1).
    """
    pop, n = K.shape
    order = np.argsort(-K, axis=1)
    X = np.zeros((pop, n), dtype=np.int8)
    for r in range(pop):
        l

…

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