Python多目标优化Skill python-multiobjective-optimization

---

name: python-multiobjective-optimization description: Python中多目标优化的专家指导 - 包括Pareto最优性、进化算法（NSGA-II、NSGA-III、MOEA/D）、标量化方法、Pareto前沿分析，以及使用pymoo、platypus和DEAP的实现 allowed-tools: [“*”]

Python多目标优化

概述

掌握多目标优化，其中必须同时优化多个冲突目标。与单目标优化找到一个最优解不同，多目标优化发现一个Pareto前沿 - 一组折衷解，改善一个目标会恶化另一个。

核心价值： 在竞争目标之间找到最优折衷（成本vs性能，风险vs回报，速度vs准确性），并实现跨Pareto前沿的知情决策。

何时使用

使用多目标优化当：

• 同时优化多个冲突目标时
• 需要在竞争目标之间进行折衷分析时
• 决策者必须从Pareto最优替代方案中选择时
• 没有偏好信息时不存在单一“最佳”解决方案时

示例：

• 投资组合优化：最大化回报，最小化风险
• 工程设计：最小化成本，最大化性能，最小化重量
• 制造业：最大化吞吐量，最小化缺陷，最小化能耗
• 药物设计：最大化疗效，最小化毒性，最大化生物利用度
• 车辆路径规划：最小化距离，最小化时间，最小化车辆数

不要使用时：

• 只有一个目标存在时（使用单目标优化）
• 目标可以用已知权重组合成单一指标时
• 一个目标远比其他目标重要时

关键概念

Pareto支配

解决方案 x 支配解决方案 y 如果：

1. x 在所有目标上不比 y 差
2. x 在至少一个目标上严格优于 y

def dominates(obj_x, obj_y, minimize=True):
    """
    检查obj_x是否支配obj_y（针对最小化）。

    如果obj_x支配obj_y则返回True。
    """
    if minimize:
        # 对于最小化：x支配y如果所有x_i <= y_i且至少一个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:
        # 对于最大化：x支配y如果所有x_i >= y_i且至少一个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

Pareto前沿（Pareto边界）

Pareto前沿 是所有非支配解决方案的集合 - 解决方案中不能改善一个目标而不恶化另一个。

属性：

• Pareto前沿上的所有点从数学角度看都同样“最优”
• 在Pareto解决方案中选择需要偏好信息
• 沿前沿移动会用一个目标交换另一个

乌托邦和纳迪尔点

• 乌托邦点：每个目标的最佳可能值（通常不可行）
• 纳迪尔点：Pareto最优解中每个目标的最差值

方法概述

标量化方法

将多目标问题转换为一系列单目标问题。

1. 加权和

• 用权重组合目标：最小化 w₁f₁(x) + w₂f₂(x)
• 优点：简单，可使用任何单目标求解器
• 缺点：无法找到非凸Pareto前沿，对缩放敏感

2. 加权度量（Lp-范数）

• 用权重最小化 ||f(x) - 理想||ₚ
• 优点：可以找到非凸前沿（p=∞）
• 缺点：需要理想点估计

3. ε-约束方法

• 最小化一个目标，约束其他：min f₁(x) s.t. f₂(x) ≤ ε
• 优点：找到整个Pareto前沿包括非凸区域
• 缺点：需要许多单目标求解，需要约束边界

4. 目标规划

• 最小化与目标目标的偏差
• 优点：当目标已知时直观
• 缺点：需要设置目标

进化算法

基于群体的元启发式，演化Pareto前沿近似。

1. NSGA-II（非支配排序遗传算法II）

• 最流行的多目标进化算法
• 快速非支配排序 + 拥挤距离
• 优点：稳健，经过良好测试，处理2-3个目标良好
• 缺点：目标多时（>3）性能下降

2. NSGA-III

• NSGA-II的扩展，用于多目标（4+）
• 基于参考点的选择
• 优点：可扩展到多目标
• 缺点：参数调优更复杂

3. MOEA/D（基于分解的多目标进化算法）

• 将问题分解为标量子问题
• 优点：高效，适合多目标
• 缺点：需要仔细设置分解

4. 其他算法

• SPEA2, PAES, GDE3, SMS-EMOA等

比较

方法	目标数	Pareto前沿类型	优点	缺点
加权和	2-3	仅凸	简单，快速	错过非凸区域
ε-约束	2-3	任何	完整覆盖	需要许多求解
NSGA-II	2-3	任何	稳健，流行	对多目标差
NSGA-III	4+	任何	处理多目标	调优复杂
MOEA/D	4+	任何	高效	分解设置

库选择

pymoo - 推荐

使用时：

• 需要现代、全面的多目标优化时
• 想要进化算法（NSGA-II、NSGA-III、MOEA/D）时
• 需要内置性能指标时
• 想要干净、面向对象的API时

优点：

• 积极维护
• 广泛的算法库
• 内置测试问题
• 性能指标（超体积等）
• 良好文档

安装：

pip install pymoo

platypus

使用时：

• 需要纯Python实现时
• 想要简单API时
• 轻量级依赖足迹时

优点：

• 纯Python（无编译）
• 简单接口
• 多种算法

安装：

pip install platypus-opt

DEAP

使用时：

• 需要通用进化计算框架时
• 想要定制灵活性时
• 已经在使用DEAP进行其他任务时

优点：

• 高度灵活
• 大社区
• 通用进化算法

缺点：

• API更复杂
• 对多目标专门化较少

安装：

pip install deap

scipy.optimize

使用时：

• 仅需要标量化方法时
• 想要最小依赖时
• 已经在使用scipy时

限制：

• 无内置进化算法
• 手动Pareto前沿构建
• 有限的多目标专门工具

实现模式

模式1：使用pymoo的NSGA-II

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

# 定义问题
class MyProblem(Problem):
    def __init__(self):
        super().__init__(
            n_var=2,           # 决策变量数
            n_obj=2,           # 目标数
            n_ieq_constr=0,    # 不等式约束数
            xl=np.array([0, 0]),    # 下界
            xu=np.array([1, 1])     # 上界
        )

    def _evaluate(self, x, out, *args, **kwargs):
        """
        评估群体x的目标。

        x: 形状为 (pop_size, n_var) 的数组
        out["F"]: 形状为 (pop_size, n_obj) 的目标
        out["G"]: 形状为 (pop_size, n_constr) 的约束 [可选]
        """
        # 目标1：最小化 (x1 - 0.5)^2 + (x2 - 0.5)^2
        f1 = (x[:, 0] - 0.5)**2 + (x[:, 1] - 0.5)**2

        # 目标2：最小化 (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])

# 创建问题实例
problem = MyProblem()

# 配置算法
algorithm = NSGA2(
    pop_size=100,           # 群体大小
    eliminate_duplicates=True
)

# 求解
res = minimize(
    problem,
    algorithm,
    ('n_gen', 200),         # 终止条件：200代
    seed=1,
    verbose=True
)

# 结果
print(f"Pareto解数：{len(res.F)}")
print(f"目标值（前5个）：
{res.F[:5]}")
print(f"决策变量（前5个）：
{res.X[:5]}")

# 可视化Pareto前沿
plot = Scatter()
plot.add(res.F, label="Pareto Front")
plot.show()

模式2：用于多目标的NSGA-III

from pymoo.algorithms.moo.nsga3 import NSGA3
from pymoo.util.ref_dirs import get_reference_directions

# 多目标问题（例如4个目标）
class ManyObjectiveProblem(Problem):
    def __init__(self):
        super().__init__(
            n_var=10,
            n_obj=4,           # 4个目标
            xl=-5,
            xu=5
        )

    def _evaluate(self, x, out, *args, **kwargs):
        # 定义4个冲突目标
        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])

# 为4个目标创建参考方向
ref_dirs = get_reference_directions("das-dennis", 4, n_partitions=12)

# NSGA-III算法
algorithm = NSGA3(
    ref_dirs=ref_dirs,
    pop_size=92  # 必须匹配参考方向
)

# 求解
problem = ManyObjectiveProblem()
res = minimize(
    problem,
    algorithm,
    ('n_gen', 300),
    seed=1,
    verbose=True
)

print(f"Pareto解数：{len(res.F)}")

模式3：带约束

class ConstrainedProblem(Problem):
    def __init__(self):
        super().__init__(
            n_var=2,
            n_obj=2,
            n_ieq_constr=2,    # 2个不等式约束
            xl=np.array([0, 0]),
            xu=np.array([5, 5])
        )

    def _evaluate(self, x, out, *args, **kwargs):
        # 目标
        f1 = x[:, 0]**2 + x[:, 1]**2
        f2 = (x[:, 0] - 3)**2 + (x[:, 1] - 3)**2

        # 不等式约束：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])

# 用NSGA-II求解（自动处理约束）
problem = ConstrainedProblem()
algorithm = NSGA2(pop_size=100)
res = minimize(problem, algorithm, ('n_gen', 200), seed=1)

# 检查约束违反
print(f"最大约束违反：{np.max(res.CV)}")  # 应该约等于0

模式4：加权和（标量化）

from scipy.optimize import minimize as scipy_minimize

def objectives(x):
    """返回目标元组"""
    f1 = x[0]**2 + x[1]**2
    f2 = (x[0] - 1)**2 + (x[1] - 1)**2
    return f1, f2

def scalarized_objective(x, weights):
    """目标的加权和"""
    f1, f2 = objectives(x)
    return weights[0] * f1 + weights[1] * f2

# 通过变化权重生成Pareto前沿
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
        })

# 转换为数组
F = np.array([[p['f1'], p['f2']] for p in pareto_front])
X = np.array([p['x'] for p in pareto_front])

print(f"找到 {len(pareto_front)} 个Pareto解")

模式5：ε-约束方法

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

# ε-约束：最小化f1，约束f2 <= epsilon
def solve_epsilon_constraint(epsilon):
    """求解给定epsilon约束在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

# 变化epsilon以追踪Pareto前沿
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"找到 {len(pareto_front)} 个Pareto解")

模式6：使用platypus

from platypus import NSGAII, Problem, Real

def evaluate(x):
    """
    返回目标列表的评估函数。

    注意：platypus期望返回列表的函数。
    """
    f1 = x[0]**2 + x[1]**2
    f2 = (x[0] - 1)**2 + (x[1] - 1)**2
    return [f1, f2]

# 定义问题
problem = Problem(2, 2)  # 2个变量，2个目标
problem.types[:] = [Real(0, 1), Real(0, 1)]  # 变量边界
problem.function = evaluate

# 用NSGA-II求解
algorithm = NSGAII(problem, population_size=100)
algorithm.run(10000)  # 10000次函数评估

# 提取结果
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"Pareto解数：{len(F)}")

模式7：使用DEAP

import random
from deap import base, creator, tools, algorithms
import numpy as np

# 设置DEAP
creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))  # 最小化两个
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):
    """返回目标元组的评估函数"""
    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)

# 初始化群体
pop = toolbox.population(n=100)

# 运行NSGA-II
algorithms.eaMuPlusLambda(
    pop, toolbox,
    mu=100,
    lambda_=100,
    cxpb=0.9,
    mutpb=0.1,
    ngen=100,
    verbose=False
)

# 提取Pareto前沿
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"Pareto解数：{len(F)}")

Pareto分析

非支配排序

def is_dominated(obj_a, obj_b, minimize=True):
    """检查obj_a是否被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):
    """
    从目标值中找到Pareto前沿。

    F: 形状为 (n_solutions, n_objectives) 的数组
    返回：非支配解的索引
    """
    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)

# 示例用法
F = np.random.rand(100, 2)  # 100个解，2个目标
pareto_idx = find_pareto_front(F, minimize=True)
pareto_F = F[pareto_idx]

print(f"Pareto前沿包含 {len(pareto_F)} 个解，共 {len(F)} 个")

超体积指标

测量Pareto前沿近似的质量。

from pymoo.indicators.hv import HV

# 超体积指标
ref_point = np.array([1.1, 1.1])  # 参考点（应支配所有解）
ind = HV(ref_point=ref_point)

# 计算超体积
hypervolume = ind(res.F)
print(f"超体积：{hypervolume:.4f}")

# 比较两个Pareto前沿
hv1 = ind(pareto_front_1)
hv2 = ind(pareto_front_2)
print(f"前沿1更好：{hv1 > hv2}")

生成距离（GD）

测量收敛到真实Pareto前沿的程度。

from pymoo.indicators.gd import GD

# 需要真实Pareto前沿（如果已知）
true_pf = np.array([...])  # 真实Pareto前沿

ind = GD(true_pf)
gd = ind(res.F)
print(f"生成距离：{gd:.4f}")  # 越低越好

逆生成距离（IGD）

测量收敛性和多样性。

from pymoo.indicators.igd import IGD

ind = IGD(true_pf)
igd = ind(res.F)
print(f"逆生成距离：{igd:.4f}")  # 越低越好

间距

测量沿Pareto前沿分布的均匀性。

def spacing(F):
    """
    计算间距指标。

    较低值表示更均匀分布。
    """
    distances = []
    for i in range(len(F)):
        # 找到到其他解的最小距离
        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"间距：{sp:.4f}")  # 越低越好

可视化

2D Pareto前沿

import matplotlib.pyplot as plt

def plot_pareto_front_2d(F, X=None, labels=None):
    """
    绘制2D Pareto前沿。

    F: 目标值，形状 (n, 2)
    X: 决策变量（可选）
    labels: 轴标签（可选）
    """
    fig, axes = plt.subplots(1, 2 if X is not None else 1, figsize=(12, 5))

    if X is None:
        axes = [axes]

    # 绘制目标空间
    axes[0].scatter(F[:, 0], F[:, 1], alpha=0.7)
    axes[0].set_xlabel(labels[0] if labels else '目标1')
    axes[0].set_ylabel(labels[1] if labels else '目标2')
    axes[0].set_title('Pareto前沿（目标空间）')
    axes[0].grid(True, alpha=0.3)

    # 绘制决策空间
    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集（决策空间）')
        axes[1].grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

# 用法
plot_pareto_front_2d(res.F, res.X, labels=['成本', '性能'])

3D Pareto前沿

def plot_pareto_front_3d(F, labels=None):
    """绘制3D Pareto前沿"""
    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 '目标1')
    ax.set_ylabel(labels[1] if labels else '目标2')
    ax.set_zlabel(labels[2] if labels else '目标3')
    ax.set_title('Pareto前沿')

    plt.show()

# 用法
plot_pareto_front_3d(res.F, labels=['成本', '时间', '质量'])

平行坐标（多目标）

from pymoo.visualization.pcp import PCP

# 对于多目标（4+）
plot = PCP()
plot.add(res.F, label="Pareto Front")
plot.show()

折衷可视化

def plot_tradeoff_curves(F, ref_idx=0):
    """
    绘制折衷曲线，显示目标沿Pareto前沿如何变化。

    ref_idx: 参考解索引
    """
    n_obj = F.shape[1]

    # 按第一个目标排序
    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('目标1')
        axes[i-1].set_ylabel(f'目标 {i+1}')
        axes[i-1].grid(True, alpha=0.3)
        axes[i-1].set_title(f'折衷：目标1 vs 目标 {i+1}')

    plt.tight_layout()
    plt.show()

plot_tradeoff_curves(res.F)

常见问题类型

投资组合优化

class PortfolioProblem(Problem):
    """
    多目标投资组合优化。

    目标：
    1. 最大化预期回报
    2. 最小化风险（方差）
    """
    def __init__(self, mu, Sigma):
        """
        mu: 预期回报 (n_assets,)
        Sigma: 协方差矩阵 (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,  # 无卖空
            xu=1.0   # 最大分配
        )

    def _evaluate(self, X, out, *args, **kwargs):
        # 将权重归一化到和为1
        W = X / X.sum(axis=1, keepdims=True)

        # 预期回报（最大化 -> 取负以最小化）
        returns = -(W @ self.mu)

        # 风险（方差）
        risk = np.array([w @ self.Sigma @ w for w in W])

        out["F"] = np.column_stack([returns, risk])

# 示例用法
np.random.seed(42)
n_assets = 10
mu = np.random.rand(n_assets) * 0.2  # 回报 0-20%
Sigma = np.random.rand(n_assets, n_assets)
Sigma = Sigma @ Sigma.T  # 使正定

problem = PortfolioProblem(mu, Sigma)
algorithm = NSGA2(pop_size=100)
res = minimize(problem, algorithm, ('n_gen', 200), seed=1)

# 绘制有效边界
plt.scatter(-res.F[:, 0], res.F[:, 1])  # 将回报取负转回正数
plt.xlabel('预期回报')
plt.ylabel('风险（方差）')
plt.title('有效边界')
plt.grid(True)
plt.show()

工程设计

class EngineeringDesign(Problem):
    """
    示例：梁设计

    目标：
    1. 最小化重量
    2. 最小化挠度

    约束：
    - 最大应力
    - 几何约束
    """
    def __init__(self):
        super().__init__(
            n_var=3,      # 宽度，高度，长度
            n_obj=2,      # 重量，挠度
            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]

        # 材料属性
        E = 200e9  # 杨氏模量 (Pa)
        rho = 7850  # 密度 (kg/m^3)
        F = 10000  # 施加力 (N)
        sigma_max = 250e6  # 最大允许应力 (Pa)

        # 目标
        weight = rho * w * h * L
        I = (w * h**3) / 12  # 面积二阶矩
        deflection = (F * L**3) / (3 * E * I)

        # 约束
        sigma = (F * L * h) / (2 * I)
        g1 = sigma - sigma_max  # 应力约束：sigma <= sigma_max
        g2 = deflection - 0.01  # 挠度约束：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解：{len(res.F)}")
print(f"最大约束违反：{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):
    """
    多目标特征选择。

    目标：
    1. 最大化准确度（最小化错误）
    2. 最小化特征数
    """
    def __init__(self, X, y):
        self.X = X
        self.y = y
        self.n_features = X.shape[1]

        super().__init__(
            n_var=self.n_features,
            n_obj=2,
            xl=0,
            xu=1,
            type_var=int  # 二进制：特征选择或不选
        )

    def _evaluate(self, X, out, *args, **kwargs):
        errors = []
        n_features = []

        for mask in X:
            # 统计选择特征
            n_sel = np.sum(mask)

            # 避免空特征集
            if n_sel == 0:
                errors.append(1.0)
                n_features.append(0)
                continue

            # 选择特征
            X_subset = self.X[:, mask.astype(bool)]

            # 用交叉验证评估
            clf = RandomForestClassifier(n_estimators=50, random_state=42)
            scores = cross_val_score(clf, X_subset, self.y, cv=3)
            error = 1 - np.mean(scores)

            errors.append(error)
            n_features.append(n_sel)

        out["F"] = np.column_stack([errors, n_features])

# 示例用法（注释以避免长运行时间）
# data = load_breast_cancer()
# problem = FeatureSelectionProblem(data.data, data.target)
# algorithm = NSGA2(pop_size=50)
# res = minimize(problem, algorithm, ('n_gen', 100), seed=1)

最佳实践

1. 归一化目标

class NormalizedProblem(Problem):
    def __init__(self, f1_scale, f2_scale):
        self.f1_scale = f1_scale
        self.f2_scale = f2_scale
        super().__init__(n_var=2, n_obj=2, xl=0, xu=1)

    def _evaluate(self, X, out, *args, **kwargs):
        # 原始目标
        f1 = X[:, 0]**2 + X[:, 1]**2
        f2 = (X[:, 0] - 1)**2 + (X[:, 1] - 1)**2

        # 归一化到相似尺度
        f1_norm = f1 / self.f1_scale
        f2_norm = f2 / self.f2_scale

        out["F"] = np.column_stack([f1_norm, f2_norm])

2. 设置合适的群体大小

经验法则：

• 2个目标：pop_size = 50-100
• 3个目标：pop_size = 100-200
• 4+个目标：pop_size = 200-500

# 自动确定群体大小
n_obj = problem.n_obj
if n_obj == 2:
    pop_size = 100
elif n_obj == 3:
    pop_size = 150
else:
    pop_size = 200 + 50 * (n_obj - 3)

algorithm = NSGA2(pop_size=pop_size)

3. 使用足够的代数

# 监控收敛
from pymoo.util.display import Display

class MyDisplay(Display):
    def _do(self, problem, evaluator, algorithm):
        super()._do(problem, evaluator, algorithm)
        # 自定义指标
        print(f"Gen {algorithm.n_gen}: HV = {ind.calc(algorithm.pop.get('F')):.4f}")

# 运行并监控
res = minimize(
    problem,
    algorithm,
    ('n_gen', 500),
    verbose=True,
    display=MyDisplay()
)

4. 多次运行以确保稳健性

# 用不同种子运行多次
results = []
hypervolumes = []

for seed in range(10):
    res = minimize(problem, algorithm, ('n_gen', 200), seed=seed, verbose=False)
    results.append(res)
    hypervolumes.append(ind(res.F))

# 选择最佳运行
best_idx = np.argmax(hypervolumes)
best_result = results[best_idx]

print(f"最佳超体积：{hypervolumes[best_idx]:.4f}")
print(f"平均超体积：{np.mean(hypervolumes):.4f} ± {np.std(hypervolumes):.4f}")

5. 验证Pareto前沿

def validate_pareto_front(F):
    """验证所有解都是非支配的"""
    n = len(F)
    for i in range(n):
        for j in range(n):
            if i != j:
                if is_dominated(F[i], F[j]):
                    print(f"警告：解 {i} 被 {j} 支配")
                    return False
    print("✓ 所有解都是非支配的")
    return True

validate_pareto_front(res.F)

决策制定

找到Pareto前沿后，决策者必须选择解。

膝点

具有最佳折衷的解（最大边际替代率）。

from pymoo.mcdm.pseudo_weights import PseudoWeights

# 找到膝点
weights = PseudoWeights()
knee_idx = weights.do(res.F)

print(f"膝点解：")
print(f"  目标：{res.F[knee_idx]}")
print(f"  变量：{res.X[knee_idx]}")

妥协规划

最小化到乌托邦点的距离。

def find_compromise_solution(F, p=2):
    """
    使用Lp-范数找到妥协解。

    p=1: 曼哈顿距离
    p=2: 欧几里得距离
    p=inf: 切比雪夫距离
    """
    # 乌托邦点（每个目标的最佳值）
    utopia = np.min(F, axis=0)

    # 纳迪尔点（每个目标的最差值）
    nadir = np.max(F, axis=0)

    # 归一化
    F_norm = (F - utopia) / (nadir - utopia + 1e-10)

    # 计算距离
    if p == np.inf:
        distances = np.max(F_norm, axis=1)
    else:
        distances = np.sum(F_norm**p, axis=1)**(1/p)

    # 找到最小距离
    best_idx = np.argmin(distances)

    return best_idx

compromise_idx = find_compromise_solution(res.F, p=2)
print(f"妥协解：{res.F[compromise_idx]}")

交互式筛选

def filter_solutions(F, X, constraints):
    """
    基于偏好筛选解。

    constraints: 字典，如 {'f1_max': 10, 'f2_max': 5}
    """
    mask = np.ones(len(F), dtype=bool)

    if 'f1_max' in constraints:
        mask &= F[:, 0] <= constraints['f1_max']
    if 'f2_max' in constraints:
        mask &= F[:, 1] <= constraints['f2_max']
    # 根据需要添加更多约束

    return F[mask], X[mask]

# 示例：筛选目标f1 <= 0.5的解
F_filtered, X_filtered = filter_solutions(
    res.F, res.X,
    {'f1_max': 0.5}
)
print(f"筛选到 {len(F_filtered)} 个解")

故障排除

问题	症状	解决方案
收敛差	Pareto前沿中有大间隙	增加代数，群体大小
分布不均匀	解聚集	使用NSGA-III或MOEA/D，增加多样性
约束违反	res.CV > 0	调整约束处理，可行初始化
运行时间长	收敛需数小时	减少群体大小，使用更快问题评估，并行化
非凸区域错过	加权和前沿中有间隙	使用ε-约束或进化算法
缩放问题	一个目标主导	归一化目标到相似范围

性能提示

并行化评估

from multiprocessing.pool import ThreadPool

# 用线程并行化
pool = ThreadPool(8)
res = minimize(
    problem,
    algorithm,
    ('n_gen', 200),
    seed=1,
    verbose=True,
    elementwise_runner=pool.starmap
)
pool.close()

从先前运行热启动

# 保存结果
np.save('pareto_front.npy', res.F)
np.save('pareto_set.npy', res.X)

# 加载并用作初始群体
prev_X = np.load('pareto_set.npy')
from pymoo.core.population import Population

pop = Population.new(X=prev_X)
algorithm.initialization = pop

# 继续优化
res = minimize(problem, algorithm, ('n_gen', 100), seed=1)

安装

# pymoo（推荐）
pip install pymoo

# platypus
pip install platypus-opt

# DEAP
pip install deap

# 可视化
pip install matplotlib

额外资源

• pymoo文档： https://pymoo.org/
• 多目标优化（Deb）： 经典教科书
• EMO会议： 进化多目标优化领先会议
• NSGA-II论文： Deb et al., “A fast and elitist multiobjective genetic algorithm: NSGA-II”

相关技能

• python-optimization - 单目标优化（scipy, pyomo, cvxpy）
• python-ase - 材料多目标优化
• pycse - 回归和参数估计