Expert guidance for multiobjective optimization in Python — Pareto optimality, evolutionary algorithms (NSGA-II, NSGA-III, MOEA/D), scalarization methods, Pareto front analysis, and implementation with pymoo, Platypus, and DEAP
Master multiobjective optimization where you must simultaneously optimize multiple conflicting objectives. Unlike single-objective optimization that finds one optimal solution, multiobjective optimization discovers a Pareto front - a set of trade-off solutions where improving one objective worsens another.
Core value: Find optimal trade-offs between competing objectives (cost vs. performance, risk vs. return, speed vs. accuracy) and enable informed decision-making across the Pareto frontier.
Use multiobjective optimization when:
Examples:
Don't use when:
Solution x dominates solution y if:
def dominates(obj_x, obj_y, minimize=True):
"""
Check if obj_x dominates obj_y (for minimization).
Returns True if obj_x dominates obj_y.
"""
if minimize:
# For minimization: x dominates y if all x_i <= y_i and at least one x_i < y_i
better_or_equal = all(x <= y for x, y in zip(obj_x, obj_y))
strictly_better = any(x < y for x, y in zip(obj_x, obj_y))
else:
# For maximization: x dominates y if all x_i >= y_i and at least one x_i > y_i
better_or_equal = all(x >= y for x, y in zip(obj_x, obj_y))
strictly_better = any(x > y for x, y in zip(obj_x, obj_y))
return better_or_equal and strictly_better
The Pareto front is the set of all non-dominated solutions - solutions where you cannot improve one objective without worsening another.
Properties:
Convert multiobjective problem to series of single-objective problems.
1. Weighted Sum
2. Weighted Metrics (Lp-norm)
3. ε-Constraint Method
4. Goal Programming
Population-based metaheuristics that evolve Pareto front approximations.
1. NSGA-II (Non-dominated Sorting Genetic Algorithm II)
2. NSGA-III
3. MOEA/D (Multiobjective Evolutionary Algorithm based on Decomposition)
4. Other Algorithms
| Method | Objectives | Pareto Front Type | Pros | Cons | |--------|-----------|-------------------|------|------| | Weighted Sum | 2-3 | Convex only | Simple, fast | Misses non-convex regions | | ε-Constraint | 2-3 | Any | Complete coverage | Many solves needed | | NSGA-II | 2-3 | Any | Robust, popular | Poor for many objectives | | NSGA-III | 4+ | Any | Handles many objectives | Complex tuning | | MOEA/D | 4+ | Any | Efficient | Decomposition setup |
Use when:
Strengths:
Installation:
pip install pymoo
Use when:
Strengths:
Installation:
pip install platypus-opt
Use when:
Strengths:
Cons:
Installation:
pip install deap
Use when:
Limitations:
import numpy as np
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.core.problem import Problem
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter
# Define problem
class MyProblem(Problem):
def __init__(self):
super().__init__(
n_var=2, # Number of decision variables
n_obj=2, # Number of objectives
n_ieq_constr=0, # Number of inequality constraints
xl=np.array([0, 0]), # Lower bounds
xu=np.array([1, 1]) # Upper bounds
)
def _evaluate(self, x, out, *args, **kwargs):
"""
Evaluate objectives for population x.
x: array of shape (pop_size, n_var)
out["F"]: objectives of shape (pop_size, n_obj)
out["G"]: constraints of shape (pop_size, n_constr) [optional]
"""
# Objective 1: minimize (x1 - 0.5)^2 + (x2 - 0.5)^2
f1 = (x[:, 0] - 0.5)**2 + (x[:, 1] - 0.5)**2
# Objective 2: minimize (x1 - 0.2)^2 + (x2 - 0.8)^2
f2 = (x[:, 0] - 0.2)**2 + (x[:, 1] - 0.8)**2
out["F"] = np.column_stack([f1, f2])
# Create problem instance
problem = MyProblem()
# Configure algorithm
algorithm = NSGA2(
pop_size=100, # Population size
eliminate_duplicates=True
)
# Solve
res = minimize(
problem,
algorithm,
('n_gen', 200), # Termination: 200 generations
seed=1,
verbose=True
)
# Results
print(f"Number of Pareto solutions: {len(res.F)}")
print(f"Objective values (first 5):\n{res.F[:5]}")
print(f"Decision variables (first 5):\n{res.X[:5]}")
# Visualize Pareto front
plot = Scatter()
plot.add(res.F, label="Pareto Front")
plot.show()
from pymoo.algorithms.moo.nsga3 import NSGA3
from pymoo.util.ref_dirs import get_reference_directions
# Problem with many objectives (e.g., 4)
class ManyObjectiveProblem(Problem):
def __init__(self):
super().__init__(
n_var=10,
n_obj=4, # 4 objectives
xl=-5,
xu=5
)
def _evaluate(self, x, out, *args, **kwargs):
# Define 4 conflicting objectives
f1 = np.sum(x[:, :5]**2, axis=1)
f2 = np.sum((x[:, :5] - 1)**2, axis=1)
f3 = np.sum(x[:, 5:]**2, axis=1)
f4 = np.sum((x[:, 5:] - 2)**2, axis=1)
out["F"] = np.column_stack([f1, f2, f3, f4])
# Create reference directions for 4 objectives
ref_dirs = get_reference_directions("das-dennis", 4, n_partitions=12)
# NSGA-III algorithm
algorithm = NSGA3(
ref_dirs=ref_dirs,
pop_size=92 # Must match reference directions
)
# Solve
problem = ManyObjectiveProblem()
res = minimize(
problem,
algorithm,
('n_gen', 300),
seed=1,
verbose=True
)
print(f"Number of Pareto solutions: {len(res.F)}")
class ConstrainedProblem(Problem):
def __init__(self):
super().__init__(
n_var=2,
n_obj=2,
n_ieq_constr=2, # 2 inequality constraints
xl=np.array([0, 0]),
xu=np.array([5, 5])
)
def _evaluate(self, x, out, *args, **kwargs):
# Objectives
f1 = x[:, 0]**2 + x[:, 1]**2
f2 = (x[:, 0] - 3)**2 + (x[:, 1] - 3)**2
# Inequality constraints: g(x) <= 0
g1 = x[:, 0] + x[:, 1] - 5 # x1 + x2 <= 5
g2 = -(x[:, 0] + 2*x[:, 1] - 6) # x1 + 2*x2 >= 6
out["F"] = np.column_stack([f1, f2])
out["G"] = np.column_stack([g1, g2])
# Solve with NSGA-II (handles constraints automatically)
problem = ConstrainedProblem()
algorithm = NSGA2(pop_size=100)
res = minimize(problem, algorithm, ('n_gen', 200), seed=1)
# Check constraint violations
print(f"Max constraint violation: {np.max(res.CV)}") # Should be ~0
from scipy.optimize import minimize as scipy_minimize
def objectives(x):
"""Return tuple of objectives"""
f1 = x[0]**2 + x[1]**2
f2 = (x[0] - 1)**2 + (x[1] - 1)**2
return f1, f2
def scalarized_objective(x, weights):
"""Weighted sum of objectives"""
f1, f2 = objectives(x)
return weights[0] * f1 + weights[1] * f2
# Generate Pareto front by varying weights
weights_list = np.linspace(0, 1, 11)
pareto_front = []
for w in weights_list:
weights = [w, 1 - w]
result = scipy_minimize(
scalarized_objective,
x0=[0.5, 0.5],
args=(weights,),
method='SLSQP',
bounds=[(0, 2), (0, 2)]
)
if result.success:
f1, f2 = objectives(result.x)
pareto_front.append({
'x': result.x,
'f1': f1,
'f2': f2,
'weights': weights
})
# Convert to arrays
F = np.array([[p['f1'], p['f2']] for p in pareto_front])
X = np.array([p['x'] for p in pareto_front])
print(f"Found {len(pareto_front)} Pareto solutions")
from scipy.optimize import minimize as scipy_minimize
def objective1(x):
return x[0]**2 + x[1]**2
def objective2(x):
return (x[0] - 1)**2 + (x[1] - 1)**2
# ε-constraint: minimize f1, subject to f2 <= epsilon
def solve_epsilon_constraint(epsilon):
"""Solve for given epsilon constraint on f2"""
constraint = {
'type': 'ineq',
'fun': lambda x: epsilon - objective2(x) # f2(x) <= epsilon
}
result = scipy_minimize(
objective1,
x0=[0.5, 0.5],
method='SLSQP',
constraints=[constraint],
bounds=[(0, 2), (0, 2)]
)
return result
# Vary epsilon to trace Pareto front
epsilon_values = np.linspace(0.01, 2.0, 20)
pareto_front = []
for eps in epsilon_values:
result = solve_epsilon_constraint(eps)
if result.success:
pareto_front.append({
'x': result.x,
'f1': objective1(result.x),
'f2': objective2(result.x),
'epsilon': eps
})
print(f"Found {len(pareto_front)} Pareto solutions")
from platypus import NSGAII, Problem, Real
def evaluate(x):
"""
Evaluation function returning list of objectives.
Note: platypus expects function that returns list.
"""
f1 = x[0]**2 + x[1]**2
f2 = (x[0] - 1)**2 + (x[1] - 1)**2
return [f1, f2]
# Define problem
problem = Problem(2, 2) # 2 variables, 2 objectives
problem.types[:] = [Real(0, 1), Real(0, 1)] # Variable bounds
problem.function = evaluate
# Solve with NSGA-II
algorithm = NSGAII(problem, population_size=100)
algorithm.run(10000) # 10000 function evaluations
# Extract results
F = np.array([[s.objectives[0], s.objectives[1]] for s in algorithm.result])
X = np.array([[s.variables[0], s.variables[1]] for s in algorithm.result])
print(f"Number of Pareto solutions: {len(F)}")
import random
from deap import base, creator, tools, algorithms
import numpy as np
# Setup DEAP
creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0)) # Minimize both
creator.create("Individual", list, fitness=creator.FitnessMin)
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, 0, 1)
toolbox.register("individual", tools.initRepeat, creator.Individual,
toolbox.attr_float, n=2)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
def evaluate(individual):
"""Evaluation function returning tuple of objectives"""
x = individual
f1 = x[0]**2 + x[1]**2
f2 = (x[0] - 1)**2 + (x[1] - 1)**2
return f1, f2
toolbox.register("evaluate", evaluate)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=0.1, indpb=0.2)
toolbox.register("select", tools.selNSGA2)
# Initialize population
pop = toolbox.population(n=100)
# Run NSGA-II
algorithms.eaMuPlusLambda(
pop, toolbox,
mu=100,
lambda_=100,
cxpb=0.9,
mutpb=0.1,
ngen=100,
verbose=False
)
# Extract Pareto front
pareto_front = tools.sortNondominated(pop, len(pop), first_front_only=True)[0]
F = np.array([ind.fitness.values for ind in pareto_front])
X = np.array([list(ind) for ind in pareto_front])
print(f"Number of Pareto solutions: {len(F)}")
def is_dominated(obj_a, obj_b, minimize=True):
"""Check if obj_a is dominated by obj_b"""
if minimize:
better = all(b <= a for a, b in zip(obj_a, obj_b))
strictly_better = any(b < a for a, b in zip(obj_a, obj_b))
else:
better = all(b >= a for a, b in zip(obj_a, obj_b))
strictly_better = any(b > a for a, b in zip(obj_a, obj_b))
return better and strictly_better
def find_pareto_front(F, minimize=True):
"""
Find Pareto front from objective values.
F: array of shape (n_solutions, n_objectives)
Returns: indices of non-dominated solutions
"""
n = len(F)
pareto_indices = []
for i in range(n):
dominated = False
for j in range(n):
if i != j and is_dominated(F[i], F[j], minimize):
dominated = True
break
if not dominated:
pareto_indices.append(i)
return np.array(pareto_indices)
# Example usage
F = np.random.rand(100, 2) # 100 solutions, 2 objectives
pareto_idx = find_pareto_front(F, minimize=True)
pareto_F = F[pareto_idx]
print(f"Pareto front contains {len(pareto_F)} solutions out of {len(F)}")
Measures quality of Pareto front approximation.
from pymoo.indicators.hv import HV
# Hypervolume indicator
ref_point = np.array([1.1, 1.1]) # Reference point (should dominate all solutions)
ind = HV(ref_point=ref_point)
# Calculate hypervolume
hypervolume = ind(res.F)
print(f"Hypervolume: {hypervolume:.4f}")
# Compare two Pareto fronts
hv1 = ind(pareto_front_1)
hv2 = ind(pareto_front_2)
print(f"Front 1 is better: {hv1 > hv2}")
Measures convergence to true Pareto front.
from pymoo.indicators.gd import GD
# Requires true Pareto front (if known)
true_pf = np.array([...]) # True Pareto front
ind = GD(true_pf)
gd = ind(res.F)
print(f"Generational Distance: {gd:.4f}") # Lower is better
Measures both convergence and diversity.
from pymoo.indicators.igd import IGD
ind = IGD(true_pf)
igd = ind(res.F)
print(f"Inverted Generational Distance: {igd:.4f}") # Lower is better
Measures uniformity of distribution along Pareto front.
def spacing(F):
"""
Calculate spacing metric.
Lower values indicate more uniform distribution.
"""
distances = []
for i in range(len(F)):
# Find minimum distance to other solutions
min_dist = np.min([np.linalg.norm(F[i] - F[j])
for j in range(len(F)) if i != j])
distances.append(min_dist)
d_bar = np.mean(distances)
spacing_metric = np.sqrt(np.sum((distances - d_bar)**2) / len(distances))
return spacing_metric
sp = spacing(res.F)
print(f"Spacing: {sp:.4f}") # Lower is better
import matplotlib.pyplot as plt
def plot_pareto_front_2d(F, X=None, labels=None):
"""
Plot 2D Pareto front.
F: Objective values, shape (n, 2)
X: Decision variables (optional)
labels: Axis labels (optional)
"""
fig, axes = plt.subplots(1, 2 if X is not None else 1, figsize=(12, 5))
if X is None:
axes = [axes]
# Plot objective space
axes[0].scatter(F[:, 0], F[:, 1], alpha=0.7)
axes[0].set_xlabel(labels[0] if labels else 'Objective 1')
axes[0].set_ylabel(labels[1] if labels else 'Objective 2')
axes[0].set_title('Pareto Front (Objective Space)')
axes[0].grid(True, alpha=0.3)
# Plot decision space
if X is not None:
axes[1].scatter(X[:, 0], X[:, 1], alpha=0.7)
axes[1].set_xlabel('x1')
axes[1].set_ylabel('x2')
axes[1].set_title('Pareto Set (Decision Space)')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Usage
plot_pareto_front_2d(res.F, res.X, labels=['Cost', 'Performance'])
def plot_pareto_front_3d(F, labels=None):
"""Plot 3D Pareto front"""
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(F[:, 0], F[:, 1], F[:, 2], alpha=0.6)
ax.set_xlabel(labels[0] if labels else 'Objective 1')
ax.set_ylabel(labels[1] if labels else 'Objective 2')
ax.set_zlabel(labels[2] if labels else 'Objective 3')
ax.set_title('Pareto Front')
plt.show()
# Usage
plot_pareto_front_3d(res.F, labels=['Cost', 'Time', 'Quality'])
from pymoo.visualization.pcp import PCP
# For many objectives (4+)
plot = PCP()
plot.add(res.F, label="Pareto Front")
plot.show()
def plot_tradeoff_curves(F, ref_idx=0):
"""
Plot trade-off curves showing how objectives change along Pareto front.
ref_idx: Reference solution index
"""
n_obj = F.shape[1]
# Sort by first objective
sorted_idx = np.argsort(F[:, 0])
F_sorted = F[sorted_idx]
fig, axes = plt.subplots(n_obj - 1, 1, figsize=(10, 4 * (n_obj - 1)))
if n_obj == 2:
axes = [axes]
for i in range(1, n_obj):
axes[i-1].plot(F_sorted[:, 0], F_sorted[:, i], 'o-', alpha=0.7)
axes[i-1].set_xlabel('Objective 1')
axes[i-1].set_ylabel(f'Objective {i+1}')
axes[i-1].grid(True, alpha=0.3)
axes[i-1].set_title(f'Trade-off: Obj 1 vs Obj {i+1}')
plt.tight_layout()
plt.show()
plot_tradeoff_curves(res.F)
class PortfolioProblem(Problem):
"""
Multiobjective portfolio optimization.
Objectives:
1. Maximize expected return
2. Minimize risk (variance)
"""
def __init__(self, mu, Sigma):
"""
mu: Expected returns (n_assets,)
Sigma: Covariance matrix (n_assets, n_assets)
"""
self.mu = mu
self.Sigma = Sigma
n = len(mu)
super().__init__(
n_var=n,
n_obj=2,
n_ieq_constr=0,
xl=0.0, # No short selling
xu=1.0 # Maximum allocation
)
def _evaluate(self, X, out, *args, **kwargs):
# Normalize weights to sum to 1
W = X / X.sum(axis=1, keepdims=True)
# Expected return (maximize -> negate for minimization)
returns = -(W @ self.mu)
# Risk (variance)
risk = np.array([w @ self.Sigma @ w for w in W])
out["F"] = np.column_stack([returns, risk])
# Example usage
np.random.seed(42)
n_assets = 10
mu = np.random.rand(n_assets) * 0.2 # Returns 0-20%
Sigma = np.random.rand(n_assets, n_assets)
Sigma = Sigma @ Sigma.T # Make positive definite
problem = PortfolioProblem(mu, Sigma)
algorithm = NSGA2(pop_size=100)
res = minimize(problem, algorithm, ('n_gen', 200), seed=1)
# Plot efficient frontier
plt.scatter(-res.F[:, 0], res.F[:, 1]) # Negate return back to positive
plt.xlabel('Expected Return')
plt.ylabel('Risk (Variance)')
plt.title('Efficient Frontier')
plt.grid(True)
plt.show()
class EngineeringDesign(Problem):
"""
Example: Beam design
Objectives:
1. Minimize weight
2. Minimize deflection
Constraints:
- Maximum stress
- Geometric constraints
"""
def __init__(self):
super().__init__(
n_var=3, # width, height, length
n_obj=2, # weight, deflection
n_ieq_constr=2,
xl=np.array([0.1, 0.1, 1.0]),
xu=np.array([1.0, 1.0, 10.0])
)
def _evaluate(self, X, out, *args, **kwargs):
w, h, L = X[:, 0], X[:, 1], X[:, 2]
# Material properties
E = 200e9 # Young's modulus (Pa)
rho = 7850 # Density (kg/m^3)
F = 10000 # Applied force (N)
sigma_max = 250e6 # Maximum allowable stress (Pa)
# Objectives
weight = rho * w * h * L
I = (w * h**3) / 12 # Second moment of area
deflection = (F * L**3) / (3 * E * I)
# Constraints
sigma = (F * L * h) / (2 * I)
g1 = sigma - sigma_max # Stress constraint: sigma <= sigma_max
g2 = deflection - 0.01 # Deflection constraint: deflection <= 0.01 m
out["F"] = np.column_stack([weight, deflection])
out["G"] = np.column_stack([g1, g2])
problem = EngineeringDesign()
algorithm = NSGA2(pop_size=100)
res = minimize(problem, algorithm, ('n_gen', 300), seed=1)
print(f"Pareto solutions found: {len(res.F)}")
print(f"Max constraint violation: {np.max(res.CV):.2e}")
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestClassifier
class FeatureSelectionProblem(Problem):
"""
Multiobjective feature selection.
Objec
<!-- Content truncated for initial SEO render. Open the source file tab for the full file. -->
Search for places (restaurants, cafes, etc.) via Google Places API proxy on localhost.
Interact with GitHub using the `gh` CLI. Use `gh issue`, `gh pr`, `gh run`, and `gh api` for issues, PRs, CI runs, and advanced queries.
Create or update AgentSkills. Use when designing, structuring, or packaging skills with scripts, references, and assets.
Start voice calls via the OpenClaw voice-call plugin.
Notion API for creating and managing pages, databases, and blocks.
Gemini CLI for one-shot Q&A, summaries, and generation.
Category:developer