Pymoo

pip install pymoo numpy matplotlib

import numpy as np
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.optimize import minimize

class SimpleBiObjective(Problem):
    def __init__(self):
        super().__init__(n_var=2, n_obj=2, xl=np.array([-2, -2]), xu=np.array([2, 2]))

    def _evaluate(self, X, out, *args, **kwargs):
        f1 = X[:, 0] ** 2 + X[:, 1] ** 2
        f2 = (X[:, 0] - 1) ** 2 + X[:, 1] ** 2
        out["F"] = np.column_stack([f1, f2])

algorithm = NSGA2(pop_size=100)
res = minimize(SimpleBiObjective(), algorithm, ("n_gen", 200), seed=1, verbose=False)
print(f"Pareto front size: {len(res.F)}")
print(f"Objective range: F1=[{res.F[:,0].min():.3f}, {res.F[:,0].max():.3f}]")

import numpy as np
from pymoo.core.problem import Problem, ElementwiseProblem

# Vectorized problem (preferred for performance)
class ZDT1(Problem):
    """ZDT1 benchmark: 30 variables, 2 objectives, known Pareto front."""
    def __init__(self):
        super().__init__(n_var=30, n_obj=2, xl=0.0, xu=1.0)

    def _evaluate(self, X, out, *args, **kwargs):
        f1 = X[:, 0]
        g = 1 + 9 * X[:, 1:].mean(axis=1)
        f2 = g * (1 - np.sqrt(f1 / g))
        out["F"] = np.column_stack([f1, f2])

# Elementwise problem with inequality constraints
class ConstrainedProblem(ElementwiseProblem):
    def __init__(self):
        super().__init__(n_var=2, n_obj=1, n_ieq_constr=2,
                         xl=np.array([-5, -5]), xu=np.array([5, 5]))

    def _evaluate(self, x, out, *args, **kwargs):
        out["F"] = (x[0] - 1) ** 2 + (x[1] - 2) ** 2  # objective
        out["G"] = np.array([
            x[0] + x[1] - 2,   # g1 <= 0
            x[0] ** 2 - x[1],  # g2 <= 0
        ])

print(f"ZDT1: {ZDT1().n_var} vars, {ZDT1().n_obj} objectives")

# Mixed-variable problem: some integer, some real
from pymoo.core.variable import Real, Integer, Choice

class MixedProblem(ElementwiseProblem):
    def __init__(self):
        vars = {
            "x": Real(bounds=(-2, 2)),
            "n": Integer(bounds=(1, 10)),
        }
        super().__init__(vars=vars, n_obj=1)

    def _evaluate(self, X, out, *args, **kwargs):
        x, n = X["x"], X["n"]
        out["F"] = (x - n) ** 2

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.algorithms.moo.nsga3 import NSGA3
from pymoo.algorithms.moo.moead import MOEAD
from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.algorithms.soo.nonconvex.de import DE
from pymoo.util.ref_dirs import get_reference_directions

# NSGA-II: best for 2-3 objectives, most widely used
nsga2 = NSGA2(pop_size=100)

# NSGA-III: designed for 3+ objectives; needs reference directions
ref_dirs = get_reference_directions("das-dennis", 3, n_partitions=12)  # ~91 dirs
nsga3 = NSGA3(pop_size=len(ref_dirs), ref_dirs=ref_dirs)

# MOEA/D: decomposition-based, good for many objectives
moead = MOEAD(ref_dirs=ref_dirs, n_neighbors=15, prob_neighbor_mating=0.7)

# GA: single-objective genetic algorithm
ga = GA(pop_size=100)

# DE: Differential Evolution, good for continuous problems
de = DE(pop_size=100, variant="DE/rand/1/bin", CR=0.9, F=0.8)

print("Algorithms initialized")

from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PM
from pymoo.operators.crossover.pntx import TwoPointCrossover
from pymoo.operators.mutation.bitflip import BitflipMutation
from pymoo.operators.sampling.rnd import FloatRandomSampling, BinaryRandomSampling

# Real-valued: Simulated Binary Crossover + Polynomial Mutation (defaults for NSGA-II)
alg_real = NSGA2(
    pop_size=100,
    sampling=FloatRandomSampling(),
    crossover=SBX(prob=0.9, eta=15),     # eta: distribution index (higher = closer to parents)
    mutation=PM(eta=20),                  # eta: higher = smaller perturbation
    eliminate_duplicates=True
)

# Binary encoding
alg_bin = GA(
    pop_size=50,
    sampling=BinaryRandomSampling(),
    crossover=TwoPointCrossover(),
    mutation=BitflipMutation(prob=0.02),
)

print("Custom operators configured")

from pymoo.termination.default import DefaultMultiObjectiveTermination
from pymoo.termination import get_termination

# Simple: fixed number of generations or evaluations
term_gen = get_termination("n_gen", 500)      # stop after 500 generations
term_eval = get_termination("n_eval", 10000)  # stop after 10,000 function evaluations

# Convergence-based (recommended for multi-objective)
term_conv = DefaultMultiObjectiveTermination(
    xtol=1e-8,      # design space tolerance
    cvtol=1e-6,     # constraint violation tolerance
    ftol=0.0025,    # objective space tolerance
    period=30,      # check every 30 generations
    n_max_gen=500,  # hard limit
    n_max_evals=100_000,
)

print("Termination criteria set")

from pymoo.optimize import minimize
import numpy as np

problem = ZDT1()
algorithm = NSGA2(pop_size=100)
res = minimize(problem, algorithm, ("n_gen", 200), seed=42, verbose=False)

# Access results
print(f"Pareto front solutions: {len(res.F)}")
print(f"Objective values (first 3):\n{res.F[:3]}")
print(f"Decision variables (first 3):\n{res.X[:3]}")
print(f"Algorithm generations: {res.algorithm.n_gen}")

# Filter for feasibility (if constraints exist)
if res.G is not None:
    feasible = (res.G <= 0).all(axis=1)
    print(f"Feasible solutions: {feasible.sum()}/{len(feasible)}")

# Performance indicators
from pymoo.indicators.hv import HV
from pymoo.indicators.igd import IGD

ref_point = np.array([1.1, 1.1])  # reference point for HV (must dominate all solutions)
hv = HV(ref_point=ref_point)
print(f"Hypervolume indicator: {hv(res.F):.4f}")

import matplotlib.pyplot as plt
from pymoo.visualization.scatter import Scatter

# Scatter plot for 2D/3D Pareto fronts
plot = Scatter(title="ZDT1 Pareto Front")
plot.add(res.F, color="blue", label="NSGA-II result")
plot.show()

# Manual matplotlib plot
fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(res.F[:, 0], res.F[:, 1], s=10, color="steelblue", alpha=0.8)
ax.set_xlabel("Objective 1 (f₁)")
ax.set_ylabel("Objective 2 (f₂)")
ax.set_title("Pareto Front — ZDT1")
plt.tight_layout()
plt.savefig("pareto_front.pdf", bbox_inches="tight")
print("Saved pareto_front.pdf")

# Parallel Coordinate Plot for 3+ objectives
from pymoo.visualization.pcp import PCP

# Generate 3-objective result for visualization
from pymoo.problems import get_problem
dtlz2 = get_problem("dtlz2")
ref_dirs = get_reference_directions("das-dennis", 3, n_partitions=12)
res3 = minimize(dtlz2, NSGA3(pop_size=len(ref_dirs), ref_dirs=ref_dirs),
                ("n_gen", 200), seed=1)

pcp = PCP(title="DTLZ2 — 3 Objectives", labels=["f1", "f2", "f3"])
pcp.add(res3.F)
pcp.show()

import numpy as np
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.optimize import minimize
import matplotlib.pyplot as plt

# Beam design: minimize weight and minimize deflection
class BeamDesign(Problem):
    """
    Variables: x[0] = width (0.1–5 cm), x[1] = height (0.5–10 cm)
    Obj 1: minimize cross-sectional area (weight proxy)
    Obj 2: minimize deflection (1/I, where I = bh³/12)
    """
    def __init__(self):
        super().__init__(n_var=2, n_obj=2,
                         xl=np.array([0.1, 0.5]),
                         xu=np.array([5.0, 10.0]))

    def _evaluate(self, X, out, *args, **kwargs):
        b, h = X[:, 0], X[:, 1]
        area = b * h                      # objective 1: area (minimize)
        I = b * h**3 / 12
        deflection = 1 / I               # objective 2: deflection (minimize)
        out["F"] = np.column_stack([area, deflection])

res = minimize(BeamDesign(), NSGA2(pop_size=100), ("n_gen", 300), seed=1)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(res.F[:, 0], res.F[:, 1], s=15, c="steelblue")
axes[0].set_xlabel("Cross-sectional area (weight)")
axes[0].set_ylabel("Deflection (1/I)")
axes[0].set_title("Pareto Front")

axes[1].scatter(res.X[:, 0], res.X[:, 1], s=15, c="coral")
axes[1].set_xlabel("Width b (cm)")
axes[1].set_ylabel("Height h (cm)")
axes[1].set_title("Design Space")

plt.tight_layout()
plt.savefig("beam_design.pdf", bbox_inches="tight")
print(f"Pareto solutions: {len(res.F)}")

import numpy as np
from pymoo.core.problem import Problem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.algorithms.moo.nsga3 import NSGA3
from pymoo.optimize import minimize
from pymoo.core.callback import Callback
from pymoo.indicators.hv import HV
from pymoo.util.ref_dirs import get_reference_directions

class HVCallback(Callback):
    def __init__(self, ref_point):
        super().__init__()
        self.hv_indicator = HV(ref_point=ref_point)
        self.history = []

    def notify(self, algorithm):
        F = algorithm.opt.get("F")
        if F is not None and len(F) > 0:
            self.history.append(self.hv_indicator(F))

problem = ZDT1()
ref_point = np.array([1.1, 1.1])

results = {}
for name, alg in [("NSGA-II", NSGA2(pop_size=100))]:
    cb = HVCallback(ref_point)
    res = minimize(problem, alg, ("n_gen", 200), callback=cb, seed=42)
    results[name] = {"res": res, "hv": cb.history}
    print(f"{name}: final HV = {cb.history[-1]:.4f}")

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 4))
for name, data in results.items():
    ax.plot(data["hv"], label=name)
ax.set_xlabel("Generation")
ax.set_ylabel("Hypervolume")
ax.legend()
plt.tight_layout()
plt.savefig("hv_convergence.pdf", bbox_inches="tight")