1. Where NSGA-II Comes From: A Brief History
Multi-objective optimization is not new. Vilfredo Pareto introduced the optimality concept that bears his name in 1896, but it took nearly 90 years for evolutionary algorithms to exploit it. The first attempt was VEGA (Vector Evaluated Genetic Algorithm) by Schaffer in 1984, which simply divided the population into subgroups, each selecting according to a different objective. It worked, but tended to converge to the extremes of the front without exploring intermediate trade-offs.
Goldberg (1989) proposed the key idea: rank solutions by Pareto dominance. Building on this, Srinivas and Deb published the original NSGA in 1994. It worked well, but the non-dominated sorting had complexity O(MN3), making it impractical for large populations. Additionally, it used fitness sharing to maintain diversity — a parameter that required manual tuning.
In 2002, Deb, Pratap, Agarwal and Meyarivan published NSGA-II in IEEE Transactions on Evolutionary Computation, solving both problems: they reduced the sorting to O(MN2) with the fast non-dominated sort algorithm, and replaced fitness sharing with crowding distance, eliminating all niche parameters. The paper has accumulated over 40,000 citations and NSGA-II remains, two decades later, the most widely used multi-objective genetic algorithm in engineering.
2. Formalizing the Multi-Objective Problem
A multi-objective optimization problem (MOP) is defined as:
min F(x) = [f1(x), f2(x), …, fM(x)]T
subject to gj(x) ≤ 0, j = 1, 2, …, J
hk(x) = 0, k = 1, 2, …, K
xiL ≤ xi ≤ xiU, i = 1, 2, …, n
Where x is the vector of decision variables in the design space Rn, F(x) is the vector of objective functions in the objective space RM, and constraints g and h define the feasible region. Unlike classical optimization, we are not looking for a single optimal point, but the Pareto set: all solutions where you cannot improve one objective without worsening at least one other.
Key concepts (assuming minimization in all objectives):
- Pareto Dominance: x1 dominates x2 if fi(x1) ≤ fi(x2) for all i AND there exists at least one j such that fj(x1) < fj(x2).
- Non-dominated Solution: No one in the population dominates it. It is a candidate to belong to the Pareto front.
- Pareto-Optimal Set: All non-dominated solutions in the entire feasible space (ideal, unreachable in practice).
- Pareto Front: Image of the Pareto-optimal set in objective space. This is what NSGA-II approximates.
Concrete example: we design a wing with 3 variables (chord, thickness, angle) and 2 objectives (minimize Cd, maximize Cl). Design A (Cd=0.012, Cl=0.80) dominates design B (Cd=0.014, Cl=0.80) because it has less drag with the same lift. But A (0.012, 0.80) does not dominate C (0.010, 0.75): A has worse Cd but better Cl, C is the opposite. Both are Pareto-optimal and belong to the front.
3. Fast Non-Dominated Sort: Ordering Without Dominating
This is the computational heart of NSGA-II. For each solution p we compute two structures:
- np: domination counter: how many solutions dominate p.
- Sp: set of solutions that p dominates.
Algorithm:
- For each pair (p, q) with p ≠ q, if p dominates q, add q to Sp. If q dominates p, increment np.
- Solutions with np = 0 form Front 1 (the non-dominated ones).
- For each solution p in the current Front, traverse Sp: decrement nq of each dominated q. If nq reaches 0, q moves to the next front.
- Repeat until all solutions are classified.
Complexity: O(MN2). The double comparison loop (step 1) costs O(MN2), and propagation across fronts is O(N2). The original NSGA (1994) did a full sort on each front, resulting in O(MN3). For N=200, the difference is 200x faster: that is why NSGA-II made large populations viable.
4. Crowding Distance: Diversity Without Parameters
The Pareto front has infinitely many solutions, but our population is finite: we need a criterion to keep the most "representative" ones. NSGA-II uses crowding distance, an estimate of the local density around each solution based on the perimeter of the hypercube formed by its nearest neighbors in each objective.
Calculation for a non-dominated front:
- Initialize distance di = 0 for all solutions in the front.
- For each objective m:
- Sort the front solutions by fm (ascending).
- Set d[first] = d[last] = ∞. This guarantees that the extremes of the front are always preserved.
- For intermediate solutions: di += (fmi+1 - fmi-1) / (fmmax - fmmin).
The normalization by (fmax - fmin) makes the metric invariant to the scale of each objective. A large distance indicates the solution is in a sparsely explored region and is valuable for the diversity of the front.
Total cost of NSGA-II per generation: O(MN2) from non-dominated sort + O(MN log N) from crowding distance = O(MN2). For N=200 and M=3, that is ~120K operations per generation. With 500 generations and each objective evaluation taking 5 minutes in CFD, the total time is dominated by the 100,000 evaluations (= 200 × 500), not by the algorithm itself.
5. The Comparison Operator: Crowded-Comparison
NSGA-II uses a special binary tournament operator, denoted ≺n:
i ≺n j if (ranki < rankj) or (ranki = rankj and di > dj)
Where rank is the Pareto front index (1 = best, 2 = next, etc.) and d is the crowding distance. In plain English: first prefer solutions from better fronts, and in case of a tie, prefer those in less crowded regions. This operator guides all selection in the algorithm (tournaments for crossover and for the next generation).
6. SBX: Simulated Binary Crossover
Simulated Binary Crossover, proposed by Deb and Agrawal in 1995, emulates the behavior of one-point crossover in binary encoding but works directly with real values. It is applied variable by variable with probability 0.5.
Given two parents x1 and x2 (with x1 < x2), the offspring are generated using a spread factor β:
c1 = 0.5 [(1+β)x1 + (1-β)x2]
c2 = 0.5 [(1-β)x1 + (1+β)x2]
Where β is sampled from a distribution controlled by the parameter ηc (typical: 15–20):
β = (2u)1/(ηc+1) if u ≤ 0.5
β = [1/(2-2u)]1/(ηc+1) if u > 0.5
Where u ~ U(0,1). High values of ηc produce offspring very close to the parents (exploitation), low values allow exploration further away. After generating offspring, they are clipped to the [xL, xU] range of each variable. The typical crossover probability is 0.9; the remaining 10% of individuals are cloned directly.
7. Polynomial Mutation
After crossover, each variable of each offspring can mutate with probability pm = 1/n (where n is the number of decision variables). If a variable x mutates:
x' = x + δ(xU - xL)
Where δ is calculated with parameter ηm (typical: 20):
δ = (2u)1/(ηm+1) - 1 if u ≤ 0.5
δ = 1 - [2(1-u)]1/(ηm+1) if u > 0.5
High ηm produces very small perturbations (local refinement), low ηm allows larger jumps. Both formulas use the same distribution transform: draw u ∈ U(0,1) and apply the inverse CDF to obtain δ with the desired polynomial tail.
8. Elitism: (μ + λ) with Pareto Fronts
NSGA-II is strongly elitist. In each generation, after generating Qt (N offspring) from Pt (N parents):
- Combine into Rt = Pt ∪ Qt (2N individuals).
- Apply fast non-dominated sort to Rt, producing fronts F1, F2, …
- Pt+1 is filled with complete fronts (F1 first, then F2, ...) as long as they fit.
- When reaching the last front that does not fit entirely, select the best solutions according to crowded-comparison ≺n.
This guarantees that the best solution found is never lost between generations. It is a (μ + λ) scheme where both parents and offspring compete to survive, unlike the (μ, λ) scheme where parents die. In problems with constraints, the comparison is modified so that any feasible solution dominates any infeasible one, and among infeasible solutions, the one with the smallest aggregate constraint violation wins.
9. Python Implementation from Scratch
Below is a complete NSGA-II implementation in pure Python, valid for any number of variables and objectives. It only requires NumPy. The code is structured as a class with independent methods for each component, making it easy to understand and extend.
import numpy as np
from dataclasses import dataclass
@dataclass
class Individual:
variables: np.ndarray
objectives: np.ndarray
rank: int = 0
crowding_distance: float = 0.0
constraints_violation: float = 0.0
class NSGA2:
def __init__(self, n_vars, n_obj, bounds, objectives_fn,
pop_size=200, max_gen=500, pc=0.9, pm=None,
eta_c=20, eta_m=20, n_constraints=0,
mutation_fn=None, crossover_fn=None):
self.n_vars = n_vars
self.n_obj = n_obj
self.bounds = np.array(bounds) # shape (n_vars, 2)
self.objectives_fn = objectives_fn
self.pop_size = pop_size
self.max_gen = max_gen
self.pc = pc
self.pm = pm if pm is not None else 1.0 / n_vars
self.eta_c = eta_c
self.eta_m = eta_m
self.n_constraints = n_constraints
self.mutation_fn = mutation_fn
self.crossover_fn = crossover_fn
self.population = []
self.fronts = []
self.history = [] # (gen, population snapshot)
def _initialize(self):
"""Initialize random population within bounds."""
self.population = []
for _ in range(self.pop_size):
vars_arr = np.random.uniform(
self.bounds[:, 0], self.bounds[:, 1]
)
obj = self._evaluate(vars_arr)
self.population.append(
Individual(variables=vars_arr, objectives=obj)
)
def _evaluate(self, x):
"""Evaluate objectives and constraints. Penalize infeasible."""
result = self.objectives_fn(x)
if self.n_constraints > 0:
obj = np.array(result[:self.n_obj])
cons = np.array(result[self.n_obj:])
else:
obj = np.array(result)
cons = np.array([])
return obj, cons
def dominates(self, ind1, ind2):
"""Does ind1 dominate ind2? (minimization in all objectives)."""
cv1 = ind1.constraints_violation
cv2 = ind2.constraints_violation
if cv1 < cv2:
return True
if cv2 < cv1:
return False
# Both feasible or same violation: compare objectives
o1, o2 = ind1.objectives, ind2.objectives
better = False
for i in range(self.n_obj):
if o1[i] > o2[i]:
return False
if o1[i] < o2[i]:
better = True
return better
def fast_non_dominated_sort(self, population):
"""Classify population into Pareto fronts."""
n = len(population)
S = [[] for _ in range(n)]
n_p = [0] * n
fronts = [[]]
for p_idx in range(n):
for q_idx in range(p_idx + 1, n):
p, q = population[p_idx], population[q_idx]
if self.dominates(p, q):
S[p_idx].append(q_idx)
n_p[q_idx] += 1
elif self.dominates(q, p):
S[q_idx].append(p_idx)
n_p[p_idx] += 1
for i in range(n):
if n_p[i] == 0:
population[i].rank = 0
fronts[0].append(i)
front_idx = 0
while front_idx < len(fronts) and fronts[front_idx]:
next_front = []
for p_idx in fronts[front_idx]:
for q_idx in S[p_idx]:
n_p[q_idx] -= 1
if n_p[q_idx] == 0:
population[q_idx].rank = front_idx + 1
next_front.append(q_idx)
front_idx += 1
if next_front:
fronts.append(next_front)
if not fronts[-1]:
fronts.pop()
return fronts
def crowding_distance(self, front_indices):
"""Compute crowding distance for a front (by index)."""
if len(front_indices) <= 2:
for idx in front_indices:
self.population[idx].crowding_distance = float('inf')
return
for idx in front_indices:
self.population[idx].crowding_distance = 0.0
for m in range(self.n_obj):
sorted_front = sorted(
front_indices,
key=lambda idx: self.population[idx].objectives[m]
)
f_min = self.population[sorted_front[0]].objectives[m]
f_max = self.population[sorted_front[-1]].objectives[m]
if f_max == f_min:
continue
self.population[sorted_front[0]].crowding_distance = float('inf')
self.population[sorted_front[-1]].crowding_distance = float('inf')
for i in range(1, len(sorted_front) - 1):
prev_obj = self.population[sorted_front[i-1]].objectives[m]
next_obj = self.population[sorted_front[i+1]].objectives[m]
self.population[sorted_front[i]].crowding_distance += \
(next_obj - prev_obj) / (f_max - f_min)
def sbx_crossover(self, parent1, parent2):
"""Simulated Binary Crossover."""
x1, x2 = parent1.variables.copy(), parent2.variables.copy()
child1, child2 = x1.copy(), x2.copy()
for i in range(self.n_vars):
if np.random.random() < 0.5:
if abs(x2[i] - x1[i]) <= 1e-14:
continue
if x1[i] < x2[i]:
lo, hi = x1[i], x2[i]
else:
lo, hi = x2[i], x1[i]
xL, xU = self.bounds[i]
u = np.random.random()
if u <= 0.5:
beta = (2 * u) ** (1.0 / (self.eta_c + 1))
else:
beta = (1.0 / (2 * (1 - u))) ** (1.0 / (self.eta_c + 1))
c1 = 0.5 * ((1 + beta) * lo + (1 - beta) * hi)
c2 = 0.5 * ((1 - beta) * lo + (1 + beta) * hi)
c1 = np.clip(c1, xL, xU)
c2 = np.clip(c2, xL, xU)
if x1[i] < x2[i]:
child1[i], child2[i] = c1, c2
else:
child2[i], child1[i] = c1, c2
return child1, child2
def polynomial_mutation(self, child):
"""Polynomial mutation variable by variable."""
mutant = child.copy()
for i in range(self.n_vars):
if np.random.random() < self.pm:
xL, xU = self.bounds[i]
u = np.random.random()
if u <= 0.5:
delta = (2 * u) ** (1.0 / (self.eta_m + 1)) - 1
else:
delta = 1 - (2 * (1 - u)) ** (1.0 / (self.eta_m + 1))
mutant[i] += delta * (xU - xL)
mutant[i] = np.clip(mutant[i], xL, xU)
return mutant
def crowded_comparison(self, ind1, ind2):
"""Crowded comparison operator."""
if ind1.rank < ind2.rank:
return -1
if ind1.rank > ind2.rank:
return 1
if ind1.crowding_distance > ind2.crowding_distance:
return -1
if ind1.crowding_distance < ind2.crowding_distance:
return 1
return 0
def tournament_selection(self, pool_indices):
"""Binary tournament selection using crowded-comparison."""
best_idx = pool_indices[0]
for idx in pool_indices[1:]:
if self.crowded_comparison(
self.population[best_idx], self.population[idx]
) > 0:
best_idx = idx
return best_idx
def evolve(self):
"""Main NSGA-II loop."""
self._initialize()
for gen in range(self.max_gen):
# Step 1: Non-dominated sort of current population
fronts = self.fast_non_dominated_sort(self.population)
self.fronts = fronts
# Step 2: Crowding distance for each front
for front in fronts:
self.crowding_distance(front)
# Step 3: Save snapshot (front 1 only)
f1_individuals = [self.population[i].objectives.copy()
for i in fronts[0]]
self.history.append((gen, f1_individuals))
# Step 4: Generate offspring Q_t (tournament + SBX + mutation)
offspring = []
pop_indices = list(range(self.pop_size))
while len(offspring) < self.pop_size:
i1 = self.tournament_selection(
np.random.choice(pop_indices, size=2, replace=False)
)
i2 = self.tournament_selection(
np.random.choice(pop_indices, size=2, replace=False)
)
if np.random.random() < self.pc:
if self.crossover_fn:
c1_vars, c2_vars = self.crossover_fn(
self.population[i1].variables,
self.population[i2].variables
)
else:
c1_vars, c2_vars = self.sbx_crossover(
self.population[i1], self.population[i2]
)
else:
c1_vars = self.population[i1].variables.copy()
c2_vars = self.population[i2].variables.copy()
if self.mutation_fn:
c1_vars = self.mutation_fn(c1_vars)
c2_vars = self.mutation_fn(c2_vars)
else:
c1_vars = self.polynomial_mutation(c1_vars)
c2_vars = self.polynomial_mutation(c2_vars)
obj1, cons1 = self._evaluate(c1_vars)
obj2, cons2 = self._evaluate(c2_vars)
ind1 = Individual(c1_vars, obj1)
ind1.constraints_violation = np.sum(np.maximum(0, cons1))
ind2 = Individual(c2_vars, obj2)
ind2.constraints_violation = np.sum(np.maximum(0, cons2))
offspring.append(ind1)
if len(offspring) < self.pop_size:
offspring.append(ind2)
# Step 5: Elitism: P_t U Q_t -> best N survive
combined = self.population + offspring
all_fronts = self.fast_non_dominated_sort(combined)
self.population = []
for front in all_fronts:
if len(self.population) + len(front) <= self.pop_size:
self.population.extend([combined[i] for i in front])
else:
remaining = self.pop_size - len(self.population)
self.crowding_distance(front)
front_sorted = sorted(
front,
key=lambda i: combined[i].crowding_distance,
reverse=True
)
self.population.extend(
[combined[i] for i in front_sorted[:remaining]]
)
break
# Return final Pareto front
final_fronts = self.fast_non_dominated_sort(self.population)
return [self.population[i] for i in final_fronts[0]]
Usage Example: ZDT1 Problem (2 objectives, 30 variables)
def zdt1(vars):
n = len(vars)
g = 1 + 9 * np.sum(vars[1:]) / (n - 1)
f1 = vars[0]
f2 = g * (1 - np.sqrt(f1 / g))
return np.array([f1, f2])
# bounds: all variables in [0, 1]
bounds = [(0.0, 1.0)] * 30
nsga = NSGA2(
n_vars=30, n_obj=2, bounds=bounds,
objectives_fn=zdt1, pop_size=200, max_gen=250
)
pareto = nsga.evolve()
print(f"Pareto front: {len(pareto)} solutions")
print(f"Best f1: {min(ind.objectives[0] for ind in pareto):.4f}")
print(f"Best f2: {min(ind.objectives[1] for ind in pareto):.4f}")
10. Constraint Handling
NSGA-II handles constraints through a simple but effective modification of the dominance operator. The rule is: a feasible solution always dominates an infeasible one. Between two infeasible solutions, the one with the smaller aggregate constraint violation dominates. Between two feasible solutions, normal Pareto dominance applies.
In the implementation, constraints_violation is computed as ∑ max(0, gj(x)) + ∑ |hk(x)|. The feasible/infeasible distinction prevents solutions with violations from being selected, but without completely discarding them from the genetic pool: they can survive if the population has not yet found the feasible region, acting as a bridge toward it.
For equality constraints hk(x) = 0, in practice they are relaxed to |hk(x)| ≤ ε with ε = 10-6 or 10-4, since exact equality is numerically impossible to satisfy with continuous variables.
11. Pareto Front Quality Metrics
Obtaining a Pareto front is not enough: you must measure its quality. The standard metrics are:
- Hypervolume (HV): volume of objective space dominated by the front, bounded by a reference point (typically the Nadir vector inflated by 10%). Measures convergence AND diversity in a single number. It is the most widely used metric in benchmarks (Zitzler et al., 2000).
- IGD (Inverted Generational Distance): average distance from points on the true Pareto front to the obtained approximation. Requires knowing the true front (only possible in synthetic benchmarks). Penalizes both missed points and false discoveries.
- Spacing (S): measures the uniformity of the solution distribution along the front. Low values indicate regular distribution. Does not measure convergence or spread.
- Coverage C(A, B): fraction of solutions in front B that are dominated by some solution in A. Useful for comparing two algorithms.
12. NSGA-III and the Future
When M ≥ 4, NSGA-II suffers from a well-known degradation: the majority of the population becomes non-dominated after a few generations and crowding distance loses discriminatory power. This is the "many-objective optimization" problem, where the very concept of dominance becomes diluted.
Deb and Jain proposed NSGA-III in 2014 as a response. It replaces crowding distance with reference points uniformly distributed over a simplex in objective space, using the Das and Dennis (1998) approach. Each solution is associated with the nearest reference point, maintaining forced diversity through niches. For M=3 with 4 divisions per axis, C(3+4-1, 4) = 15 reference points are generated; for M=10 with 5 divisions, C(14,5) = 2,002 points cover the 10-dimensional space uniformly.
Other modern extensions include MOEA/D (Zhang and Li, 2007), which decomposes the MOP into N scalar subproblems using weights and solves them simultaneously with neighbor collaboration, and SMS-EMOA (Beume et al., 2007), which uses hypervolume directly as a selection criterion. In practice, for M ≤ 3, NSGA-II remains the de facto standard.
13. OpenFOAM Integration: Optimization with CFD
One of the most powerful use cases of NSGA-II in engineering is aerodynamic optimization coupled with CFD. OpenFOAM, the most widely used open-source solver in industry, integrates with NSGA-II through an external loop where each objective function evaluation triggers a complete simulation.
Typical architecture for wing or 2D airfoil optimization:
- Design variables: n geometric parameters (e.g., 5 control points of a B-spline defining the airfoil shape, plus angle of attack).
- Geometry generation: a Python script generates the STL file or volume mesh with blockMesh/snappyHexMesh from the variables.
- CFD simulation: simpleFoam in turbulent regime (RANS with k-ω SST model). Fixed boundary conditions for all geometries.
- Post-processing: extract Cl and Cd from force files (forceCoeffs).
- Objectives: minimize Cd, maximize Cl (or minimize -Cl). Optional constraint: minimum airfoil thickness for structural integrity.
Runtime is dominated by CFD simulations. With NSGA-II, population N=100 and 200 generations, 20,000 evaluations are needed. At 5 minutes per simulation, total time is ~70 days on a single core. This is why acceleration strategies are used in practice:
- Surrogate models (Kriging, GP): train a metamodel with the first ~200 simulations and use it to pre-filter candidates.
- Parallel evaluation: NSGA-II evaluates the entire offspring population simultaneously; with a 100-node cluster, 20,000 simulations complete in ~17 hours.
- Early stopping criteria: stop CFD when residuals drop below 10-3 instead of 10-6 for the exploratory phase, refining only the promising solutions from the final front.
14. Recommended Parameters and Calibration
| Parameter | Symbol | Typical Range | Recommendation | Notes |
|---|---|---|---|---|
| Population | N | 40–500 | 200 | Larger N gives better front coverage but more evaluations |
| Generations | G | 100–2000 | 500 | Monitor hypervolume for convergence |
| Crossover prob. | pc | 0.6–1.0 | 0.9 | Values < 0.8 only in highly multimodal problems |
| Mutation prob. | pm | 1/n_vars | 1/n | The theoretical optimum according to Deb |
| SBX index | ηc | 5–50 | 20 | Higher = offspring closer to parents |
| Mutation index | ηm | 5–50 | 20 | Higher = smaller mutations |
To calibrate ηc and ηm, the practical rule is to measure the percentage of feasible offspring after crossover/mutation. If fewer than ~70% of offspring are feasible, increase ηc and reduce pm to be less aggressive with perturbations.
15. Practical Final Recommendations
- Variable and objective scaling: always normalize objectives before computing crowding distance (the code above does this implicitly with range normalization). Variables in [0,1] improve SBX numerical stability.
- Visualization in 2D/3D: use pair plots for M>3. The Pareto front in 4+ dimensions is hard to visualize; resort to 2D projections by objective pairs and parallel coordinate plots.
- Post-optimization decision making: the Pareto front is the starting point, not the end. Methods like TOPSIS or minimum distance to the ideal point (utopia point) help select a compromise solution. In practice, present 3-5 representative solutions to the engineer: best balance, best in objective 1, best in objective 2.
- Seeds: NSGA-II is stochastic. Repeat the optimization with 5-10 different seeds and report the best front by hypervolume or the envelope of all fronts combined (non-dominated sort over the union of all solutions found).
- Do not reinvent the wheel: pymoo (Blank and Deb, 2020) offers implementations of NSGA-II, NSGA-III, MOEA/D and a dozen more algorithms with included benchmarks. Use pymoo for production and the code in this article to understand the fundamentals.
Bring multi-objective optimization to your designs
At MarteriaTech Labs we integrate NSGA-II with CFD simulation, structural analysis, and parametric design to solve real engineering problems. From aerodynamic wing optimization to material selection with conflicting criteria.
Talk to us