pymoo

pymoo

概述

pymoo 通过基于种群的进化算法，为多目标优化提供统一的 API。用户通过继承 Problem 或 ElementwiseProblem 定义问题，并指定目标数（n_obj）、决策变量数（n_var）和可选的约束数（n_ieq_constr）。NSGA-II 和 NSGA-III 等算法会返回一个 Result 对象，其中包含帕累托最优种群、目标值和决策变量值。pymoo 将问题定义、算法配置、算子选择和分析分离开来，每个组件都可独立替换。

适用场景

• 优化具有两个或更多冲突目标的设计（例如，在提升性能的同时降低成本）
• 在没有梯度时，将进化算法（GA、DE、PSO）用作黑盒优化器
• 对准确率和推理时间相互制约的 ML 模型进行多目标超参数搜索
• 为投资组合优化或多准则决策分析计算帕累托前沿
• 针对特定领域的解编码（二进制、排列、实数值）定制交叉／变异算子
• 在标准测试问题（ZDT、DTLZ、CTP）上对优化算法进行基准测试
• 对于单目标、可获得梯度且平滑的优化，请改用 scipy.optimize

前置条件

• Python 包：pymoo、numpy、matplotlib
• 数据要求：目标函数和可选的约束函数；变量边界
• 环境：对于大多数问题，CPU 已足够；pymoo 核心不使用 GPU

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}]")

核心 API

模块 1：问题定义

通过继承定义优化问题。向量化评估使用 Problem（速度更快），标量评估使用 ElementwiseProblem（编写更简单）。

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

模块 2：算法选择

pymoo 提供 20 多种算法。按问题类型选择的主要算法如下：

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")

模块 3：算子（交叉与变异）

算子决定解如何演化。替换默认算子以匹配变量类型。

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")

模块 4：终止条件

控制算法何时停止。

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")

模块 5：结果分析与帕累托前沿

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}")

模块 6：可视化

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()

核心概念

帕累托支配

如果解 a 在所有目标上都不劣于 b，且至少在一个目标上严格优于 b，则 a 支配 b。帕累托前沿是非支配解的集合——不存在唯一的“最佳”解，只有取舍。NSGA-II 使用非支配排序和拥挤距离来维持具有多样性的帕累托近似解。

约束处理

pymoo 采用约束违反方法：不可行解会受到惩罚，但仍保留在种群中。无论目标值如何，违反约束 G[i] > 0 的解都会被任何可行解支配。这意味着算法会先推动种群满足可行性，再优化目标。

常见工作流

工作流 1：双目标工程设计

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)}")

工作流 2：使用回调比较算法

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")

关键参数

参数	模块/类	默认值	范围/选项	作用
pop_size	所有算法	100	50–500	每代的种群规模；越大，多样性越好，但速度越慢
n_gen	终止条件	—	100–5000	运行的最大代数
eta（交叉）	SBX	15	5–30	分布指数；越高，子代与父代越接近
eta（变异）	PM	20	5–50	扰动强度；越高，变异步长越小
CR	DE	0.9	0–1	差分进化中的交叉概率
F	DE	0.8	0.1–2.0	差分进化变异的缩放因子
n_partitions	get_reference_directions	12	4–20	NSGA-III/MOEA/D 参考方向的密度
n_neighbors	MOEAD	15	5–30	MOEA/D 权重向量选择的邻域大小
prob	SBX	0.9	0.5–1.0	对一对个体执行交叉的概率

最佳实践

1. 先分析目标函数的性能：pymoo 会调用目标函数 pop_size × n_gen 次。如果单次评估耗时超过 1 ms，请使用 pymoo.core.problem.StarmapParallelization 或 dask 并行化。先分析性能，再进行优化。

2. 2 个目标时从 NSGA-II 开始，3 个及以上目标时使用 NSGA-III：NSGA-II 是双目标问题的事实标准。对于 3 个及以上目标，拥挤距离的效果会下降，应使用 Das-Dennis 参考方向的 NSGA-III 或 MOEA/D。

3. 根据收敛情况设置终止条件，而非固定代数：DefaultMultiObjectiveTermination 会自动检测停滞。如果 Pareto 前沿提前收敛，固定的 n_gen 会浪费计算资源；如果问题较难，则可能过早终止。

4. 使用向量化的 Problem，而非 ElementwiseProblem，以提升速度：ElementwiseProblem 一次评估一个解；Problem 则通过一次 NumPy 调用评估整个种群。对于兼容 NumPy 的函数，速度可提升 10–100 倍。

5. 计算指标前先对目标进行归一化：超体积和 IGD 对目标的尺度很敏感。如果目标的单位不同（例如，质量以 kg 为单位，挠度以 m⁻¹ 为单位），请先将其归一化到 [0, 1] 再进行比较。

常见方案

方案：并行执行耗时较长的目标评估

from multiprocessing.pool import Pool
from pymoo.core.problem import StarmapParallelization

# Wrap a slow objective function with parallel evaluation
def expensive_objective(x):
    import time; time.sleep(0.01)  # simulate slow call
    return [x[0]**2 + x[1]**2, (x[0]-1)**2 + x[1]**2]

n_workers = 4
pool = Pool(n_workers)
runner = StarmapParallelization(pool.starmap)

class ParallelProblem(Problem):
    def __init__(self, runner):
        super().__init__(n_var=2, n_obj=2, xl=-2, xu=2, elementwise=True,
                         elementwise_runner=runner)
    def _evaluate(self, x, out, *args, **kwargs):
        out["F"] = expensive_objective(x)

res = minimize(ParallelProblem(runner), NSGA2(pop_size=50), ("n_gen", 50))
pool.close()
print(f"Pareto solutions: {len(res.F)}")

方案：从之前的种群重新开始

import numpy as np
from pymoo.core.population import Population

# Warm-start: use previous result's population as initial population
# Run initial optimization
res1 = minimize(ZDT1(), NSGA2(pop_size=100), ("n_gen", 100), seed=1)

# Continue from checkpoint
initial_pop = Population.new("X", res1.X)

from pymoo.algorithms.moo.nsga2 import NSGA2
alg_warmstart = NSGA2(pop_size=100, sampling=initial_pop)
res2 = minimize(ZDT1(), alg_warmstart, ("n_gen", 100), seed=1)
print(f"Continued optimization: {len(res2.F)} Pareto solutions")

方案：使用自定义适应度的单目标 GA

from pymoo.algorithms.soo.nonconvex.ga import GA
from pymoo.core.problem import ElementwiseProblem
from pymoo.optimize import minimize

class Sphere(ElementwiseProblem):
    def __init__(self):
        super().__init__(n_var=5, n_obj=1, xl=-5.0, xu=5.0)

    def _evaluate(self, x, out, *args, **kwargs):
        out["F"] = sum(xi**2 for xi in x)

res = minimize(Sphere(), GA(pop_size=100), ("n_gen", 200), seed=42)
print(f"Best solution: f={res.F[0][0]:.6f}")
print(f"Best x: {res.X}")

故障排查

问题	原因	解决方案
Pareto 前沿只有 1–2 个解	pop_size 太小或 n_gen 太少	增大 pop_size（≥100）和 n_gen（2 个目标时 ≥200）
许多代之后所有解仍不可行	约束过严，或初始采样未覆盖可行区域	通过 sampling 参数添加一个可行的种子解；为热启动放宽约束
目标值中出现 NaN 或 inf	问题定义存在数值不稳定性	在 _evaluate 中添加边界检查；除法前使用 np.clip
StarmapParallelization 挂起	Pool 未关闭；lambda/闭包不可序列化	使用 pool.close() + pool.join()；将目标函数定义为模块级函数
NSGA-III 在 2 个目标上的表现不如 NSGA-II	NSGA-III 专为 3 个及以上目标设计；在二维情况下选择压力较小	2 个目标时使用 NSGA-II；3 个及以上目标时使用 NSGA-III
收敛过早停滞	种群已收敛，缺乏多样性	增大变异的 eta（加强扰动）；增大 pop_size；改用 DE
不同运行之间的结果差异很大	随机性较高；未设置种子	在 minimize() 中设置 seed=42，以确保结果可复现

相关技能

• scipy-optimization — 针对平滑问题的单目标、基于梯度的优化
• pymatgen — 使用 pymoo 优化材料属性，进行多目标晶体结构搜索
• scikit-learn-machine-learning — 超参数调优（多目标 HPO 可使用 pymoo：准确率与延迟之间的权衡）

参考资料

• pymoo 文档 — 算法目录、API 参考和教程
• pymoo 论文：Blank 和 Deb（2020），IEEE Access — 算法说明和基准测试
• NSGA-II 论文：Deb 等（2002），IEEE TEC — 多目标 EA 奠基性论文
• NSGA-III 论文：Deb 和 Jain（2014），IEEE TEC — 面向多目标优化的扩展
• pymoo GitHub — 源代码、示例和问题跟踪器