pymoo

pymoo

概要

pymoo は、集団ベースの進化アルゴリズムを用いた多目的最適化のための統一 API を提供します。ユーザーは Problem または ElementwiseProblem をサブクラス化し、目的関数（n_obj）、設計変数（n_var）、および任意の制約（n_ieq_constr）を指定して問題を定義します。NSGA-II や NSGA-III などのアルゴリズムは、パレート最適な集団、目的関数値、設計変数値を含む Result オブジェクトを返します。pymoo では、問題定義、アルゴリズム設定、オペレーターの選択、分析を分離しており、各コンポーネントを個別に置き換えられます。

使用する場面

• 互いに競合する目的が 2 つ以上ある設計を最適化する場合（例：コストを最小化しながら性能を最大化する）
• 勾配を利用できない場合に、進化アルゴリズム（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 より悪くなく、少なくとも 1 つの目的で b より優れている場合です。パレートフロントは、他の解に支配されない解の集合です。単一の「最良」解はなく、トレードオフだけが存在します。NSGA-II は非支配ソートと混雑距離を使用して、多様なパレート近似を維持します。

制約の処理

pymoo は制約違反のアプローチを使用します。実行不可能な解にはペナルティが課されますが、集団内には残されます。制約違反 G[i] > 0 の解は、目的関数値にかかわらず、実行可能な解に支配されます。つまり、アルゴリズムはまず集団を実行可能な状態へ導き、その後、目的関数を最適化します。

一般的なワークフロー

ワークフロー 1：2 目的のエンジニアリング設計

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 回の評価に 1 ms 以上かかる場合は、pymoo.core.problem.StarmapParallelization または dask を使用して並列化してください。最適化の前にプロファイルを取得します。

2. 目的関数が 2 つの場合は NSGA-II、3 つ以上の場合は NSGA-III から始めます：NSGA-II は 2 目的問題の事実上の標準です。目的関数が 3 つ以上の場合、混雑距離の性能が低下するため、Das-Dennis 参照方向を使用する NSGA-III または MOEA/D を使用してください。

3. 固定世代数ではなく、収束状況に基づいて終了条件を設定します：DefaultMultiObjectiveTermination は停滞を自動的に検出します。パレートフロントが早期に収束した場合、固定の n_gen は計算リソースを無駄にします。また、問題が難しい場合は早く終了しすぎることがあります。

4. 速度を優先する場合は ElementwiseProblem ではなく、ベクトル化された Problem を使用します：ElementwiseProblem は一度に 1 つの解を評価しますが、Problem は NumPy の呼び出し 1 回で集団全体を評価します。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}")

トラブルシューティング

問題	原因	解決策
パレートフロントの解が 1～2 個しかない	pop_size が小さすぎるか、n_gen が少なすぎます	pop_size（≥100）と n_gen（目的関数が 2 つの場合は ≥200）を増やしてください
多くの世代を経ても、すべての解が実行可能ではない	制約が厳しすぎるか、初期サンプリングで実行可能領域を見逃しています	sampling パラメーターで実行可能なシード解を追加してください。ウォームスタートでは制約を緩和してください
目的関数の値に NaN または inf が含まれる	問題定義における数値の不安定性	_evaluate に境界チェックを追加し、除算の前に np.clip を使用してください
StarmapParallelization が停止する	Pool が閉じられていない、または lambda／クロージャーが pickle 化できません	pool.close() と pool.join() を使用してください。目的関数はモジュールレベルの関数として定義してください
目的関数が 2 つの場合、NSGA-III の性能が NSGA-II を下回る	NSGA-III は目的関数が 3 つ以上の場合に向けて設計されており、2 次元では選択圧が弱くなります	目的関数が 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 et al.（2002）、IEEE TEC — 多目的 EA の基礎となる論文
• NSGA-III の論文：Deb & Jain（2014）、IEEE TEC — 多目的拡張
• pymoo の GitHub — ソース、サンプル、イシュートラッカー