Pymoo：优化算法收敛性的实例分析

• 😈 Github 仓库： Pymoo: 基于 Python 的多目标优化仓库 (MOO, Multi-objective Optimization)。
• 🏠 Pymoo文档首页：Multi-objective Optimization in Python

✨ 注意：本文使用Pymoo版本为0.6.0。

🦄 优化算法的收敛性分析的必要性：
优化算法的收敛性是评价算法性能的一个重要指标。算法的收敛性可以为我们设计算法的终止准则(Termination Criterion)提供一定参考，收敛性分析需要通常考虑如下两个方面：

• Pareto前沿未知；
• Pareto前沿已经分析推导出来，或存在合理的近似值。

本文我们将使用IGD与Hypervolume两个性能指标分别对已知Pareto前沿与未知Pareto前沿两种情况的算法性能。

一、优化问题的定义与优化算法的结果分析

1.1 优化问题的定义与必要的算法中间信息的获取

• 1、优化问题的定义

为了使用Pymoo进行算法收敛性分析之前，我们首先应该得到一个优化问题及其分析结果。这里以文献¹中Pymoo定义的具有两个约束及双目标优化问题MyProblem及其分析结果为例，进行本文的优化算法收敛性分析过程。

另外，根据文献²这一篇文章，我们得知由于MyProblem问题相对来说比较简单，我们可以得到其Pareto前沿。

• MyProblem的Pymoo实现代码如下所示：

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

class MyProblem(ElementwiseProblem):

    def __init__(self):
        super().__init__(
            n_var        = 2,       # 两个自变量
            n_obj        = 2,       # 两个目标函数
            n_ieq_constr = 2,       # 两个不等式约束
            xl           = np.array([-2, -2]),  # 两个自变量下限
            xu           = np.array([2, 2])     # 两个自变量上限
        )

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

        g1 = 2 * (x[0] - 0.1) * (x[0] - 0.9) / 0.18
        g2 = -20 * (x[0] - 0.4) * (x[0] - 0.6) / 4.8

        out["F"] = [f1, f2]     # 计算得到的目标存放的位置
        out["G"] = [g1, g2]     # 计算得到的约束存放的位置
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

• MyProblem使用的NSGA2算法Pymoo algorithm对象实现代码如下所示：

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.operators.crossover.sbx import SBX
from pymoo.operators.mutation.pm import PM
from pymoo.operators.sampling.rnd import FloatRandomSampling

algorithm = NSGA2(
    pop_size             = 40,      # 种群规模
    n_offspring          = 10,      # 生成的子代规模
    sampling             = FloatRandomSampling(),   # 随机抽样算子
    crossover            = SBX(prob=0.9, eta=15),   # 模拟二元交叉算子
    mutation             = PM(eta=20),              # 多项式突变算子
    eliminate_duplicates = True     # 去除重复的个体
)
1
2
3
4
5
6
7
8
9
10
11
12
13

在上面定义的MyProblem及其分析算法algorithm的基础上，为了进一步检查算法分析出的最优解与理论最优解的接近程度，我们需要将计算最优解与Pareto前沿的接近程度，以便分析算法的收敛性。为了直观地感受这个过程，下面定义一个继承于MyProblem问题的MyTestProblem对象：

from pymoo.util.misc import stack

class MyTestProblem(MyProblem):

    # 计算优化问题的Pareto前沿（目标空间）
    def _calc_pareto_front(self, flatten=True, *args, **kwargs):
        # 理论上分析一个标优化问题的Pareto前沿，其实就是寻找目标函数之间的函数关系
        # 具体为什么 f1 与 f2 的关系为如下表达式，可以参考：
        # https://blog.csdn.net/weixin_37926734/article/details/127392027?spm=1001.2014.3001.5501
        f2 = lambda f1: ((f1 / 100) ** 0.5 - 1)**2
        # f1 满足Pareto集的区间
        F1_a, F1_b = np.linspace(1, 16, 300), np.linspace(36, 81, 300)
        # 根据 f1 的区间，得到对应的 f2 值
        F2_a, F2_b = f2(F1_a), f2(F1_b)

        # 由于Pareto前沿是分段函数，所以需要分开处理
        pf_a = np.column_stack([F1_a, F2_a])
        pf_b = np.column_stack([F1_b, F2_b])

        return stack(pf_a, pf_b, flatten=flatten)

    # 计算优化问题的Pareto集 (变量空间)
    def _calc_pareto_set(self, *args, **kwargs):
        # 定义优化问题的Pareto集合的本质是确定满足问题最优解的自变量区间
        # 具体的区间确定方法同样可以参考：https://blog.csdn.net/weixin_37926734/article/details/127392027?spm=1001.2014.3001.5501
        # 对于本文的`MyProblem`定义的问题其Pareto集合可以从理论分析上得到：
        # x2 = 0, x 属于 [0.1, 0.4]与[0.6, 0.9] 两个区间，因此代码如下所示：
        x1_a = np.linspace(0.1, 0.4, 50)
        x1_b = np.linspace(0.6, 0.9, 50)
        x2 = np.zeros(50)

        a, b = np.column_stack([x1_a, x2]), np.column_stack([x1_b, x2])
        return stack(a, b, flatten=flatten)

problem = MyTestProblem()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

在Pymoo中进行优化算法收敛性分析时，我们需要保存历史数据。一种方式就是通过设置save_history=True参数，即对算法对象每次迭代的结果进行深拷贝，并存入到Result对象中，对MyTestProblem进行最优化迭代，并将最优解及最优值存储到X与F中，其代码如下所示：

from pymoo.optimize import minimize

res = minimize(
    problem,
    algorithm,
    ('n_gen', 40),
    save_history=True,
    verbose=False
)

X, F = res.opt.get("X", "F")
hist = res.history
1
2
3
4
5
6
7
8
9
10
11
12

• 2、必要的算法中间信息的获取

同时，我们还需要考虑算法的中间步骤，通过如下两步实现：

• 一个Callback类用于存储算法每次迭代必要的信息。
• 使能save_history参数，使得最小化方法能够对每次迭代的算法对象进行深度拷贝。

至此，我们就可以通过下面的代码提取评价算法收敛性的必要信息的代码：

n_evals = []        # 函数迭代的次数（Corresponding number of function evaluations）
hist_F = []         # 每一代目标空间中的目标值
hist_cv = []        # 在每一代的约束冲突 (Constraint Violation)
hist_cv_avg = []    # 整个种群的平均约束冲突 (Average Constraint Violation)

for algo in hist:
    # Store the number of function evaluations
    # 存储功函数计算的次数
    n_evals.append(algo.evaluator.n_eval)

    # retrieve the optimum from the algorithm
    # 从算法中检索最优值
    opt = algo.opt

    # store the least contraint violation and the average in each population
    # 存储最新的约束冲突与种群的平均约束冲突
    hist_cv.append(opt.get("CV").min())
    hist_cv_avg.append(algo.pop.get("CV").mean())

    # filter out only the feasible and append and objective space values
    # 只过滤出可行的、附加的和目标的空间值
    feas = np.where(opt.get("feasible"))[0]
    hist_F.append(opt.get("F")[feas])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

1.2 Constraint Satisfaction (约束满足)

首先，快速获取算法是在何时得到第一个可行解(Feasible Solution)：

k = np.where(np.array(hist_cv) <= 0.0)[0].min()
print(f"At least one feasible solution in Generation {k} after {n_evals[k]} evaluations.")
1
2

代码执行结果如下图所示：

由于本文的问题不是很复杂，所以上面的代码能够很快寻找到一个最优解。然而，对有实际中的复杂优化问题，我们通常还需要执行如下所示的代码分析：

from matplotlib import pyplot as plt
# 如果最新的可行解而不是总体的平均值，可以替换为 hist_cv

vals = hist_cv_avg

k = np.where(np.array(vals) <= 0.0)[0].min()
print(f"Whole population feasible in Generation {k} after {n_evals[k]} eavaluations.")

plt.figure(figsize=(7, 5))
plt.plot(n_evals, vals, color="black", lw=0.7, label="Avg. CV of Pop")
plt.scatter(n_evals, vals, facecolor="none", edgecolor='black', marker="p")

# 绘制一条垂直于 x 轴的一条经过第一次出现可行解处
plt.axvline(n_evals[k], color="red", label="All feasible", linestyle="--")

plt.title("Convergence")
plt.xlabel("Function Evaluations")
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

代码执行结果如下图所示：

二、Pareto前沿已知的优化问题的收敛性分析

2.1 获取优化问题Pareto前沿的方法

文本中的优化问题由于相对简单，所以其Pareto前沿可以手动实现，也可以直接在问题定义中实现，以便动态分析运行，本文使用算法历史数据作为后处理的数据来源。

• 实际优化问题的Pareto前沿近似方法：

然而，对于实际问题，优化问题的Pareto前沿通常需要使用近似值。即通过多次运行算法并从所有解集中提取非支配解，以便获得Pareto前沿的近似值。如果只有一次运行，另一种方法是使用获得的非支配解集作为近似值。然而，结果只表明算法收敛到最终解集的进度。

在Pymoo的框架内部，本文定义的MyTestProblem类的实例对象problem的Pareto前沿可以通过下面的代码计算得到，并对其进行可视化处理：

# 获取Pareto前沿
# pf = problem.pareto_front(use_cache=False, flatten=True)
pf_a, pf_b = problem.pareto_front(use_cache=False, flatten=False)

# 绘制结果
plt.figure(figsize=(7, 5))

# 首先绘制优化算法计算得到的最优解
plt.scatter(F[:, 0], F[:, 1], s=30, facecolor="none", edgecolors="b", label="Solutions")

# 然后绘制理论分析得到的Pareto前沿 pf_a, pf_b
plt.plot(pf_a[:, 0], pf_a[:, 1], alpha=0.5, linewidth=2.0, color="red", label="Pareto Front")
plt.plot(pf_b[:, 0], pf_b[:, 1], alpha=0.5, linewidth=2.0, color="red")

plt.title("Objective Space")
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

代码执行结果如下所示：

2.2 Pareto前沿已知的优化算法收敛性分析

对于已知Pareto前沿的优化问题，Pymoo中可用的算法性能评价指标包括IGD/GD/IDG+/GD+(具体可以参考文献³)，本文使用IGD与IGD+两个性能指标对本文的优化问题的收敛性进行分析。

🚀 IGD优化算法收敛性评价指标说明：
（1）IGD³ (Inverted Generational Distance) 算法评价指标计算目标空间(Objective Space)中从Pareto前沿中的每个解到非支配集合中距离自己最近解的平均欧几里德距离(Average Euclidean Distance)。
（2）IGD+ (Inverted Generational Distance Plus)算法评价指标使用一个考虑支配关系的距离度量代替IGD中使用的欧几里德距离，进而实现了弱Pareto的相容性。

• 1、基于IGD性能指标的收敛性分析

from pymoo.indicators.igd import IGD

# 获取求解问题的Pareto前沿
pf = problem.pareto_front(use_cache=False, flatten=True)

# 算法性能评价指标IGD对象的定义
metric = IGD(pf, zero_to_one=True)

# 计算优化算法计算的最优解与Pareto前沿的IGD指标
igd = [metric.do(_F) for _F in hist_F]

# 绘制收敛分析结果
# 首先绘制igd性能指标的结果
plt.plot(n_evals, igd, color='black', lw=0.7, label="Avg. CV of Pop")
plt.scatter(n_evals, igd, facecolor="none", edgecolor="black", marker="p")

# 绘制一条满足收敛条件（这里假设为10^-2）的一条垂直于y轴的直线
plt.axhline(10**(-2), color="red", label="$10^{-2}$", linestyle="--")

plt.title("Convergence")
plt.xlabel("Function Evaluations")
plt.ylabel("IGD")
# y 轴使用log刻度
plt.yscale("log")
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

代码执行结果如下图所示：

• 2、基于IGD+性能指标的收敛性分析

from pymoo.indicators.igd_plus import IGDPlus

# 算法性能评价指标IGD+对象的定义
metric = IGDPlus(pf, zero_to_one=True)

# 计算优化算法计算的最优解与Pareto前沿的IGD+指标
igd_plus = [metric.do(_F) for _F in hist_F]

# 绘制收敛分析结果
# 首先绘制igd+性能指标的结果
plt.plot(n_evals, igd_plus, color="black", lw=0.7, label="Avg. CV of Pop")
plt.scatter(n_evals, igd_plus, facecolor="none", edgecolor="black", marker="p")

# 绘制一条满足收敛条件（这里假设为10^-2）的一条垂直于y轴的直线
plt.axhline(10**-2, color="red", label="$10^{-2}$", linestyle="--")

plt.title("Convergence")
plt.xlabel("Function Evaluations")
plt.ylabel("IGD+")
plt.yscale("log")
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

代码执行结果如下图所示：

三、Pareto前沿未知的优化问题的收敛性分析

如果Pareto前沿未知，我们无法知道算法是否收敛到真正的最优值，至少没有任何进一步的信息。然而，我们可以判断算法在迭代的过程中算法优化的程度，如果随着迭代的进行几乎没有任何改进的时候，则终止迭代。

🚀 此外，Pymoo还提供了针对Pareto前沿未知的优化问题收敛性分析的两种实现方法：

• （1）Hypervolume (HV)
• （2）Running Metric

3.1 基于Hypervolume的收敛性分析

Hypervolume是多目标问题优化的一个非常著名的性能指标，它符合Pareto标准，基于预定义的参考点与计算得到的最优解之间的体积。因此，Hypervolume需要定义一个参考点ref_point，并且它应大于Pareto前沿的最大值，代码如下所示：

# 沿着F的行方向，获得计算得到的目标值中的最小值与最大值
approx_idel = F.min(axis=0)
approx_nadir = F.max(axis=0)

from pymoo.indicators.hv import Hypervolume

# 算法性能评价指标Hypervolume对象的定义
metric = Hypervolume(
    ref_point=np.array([1.1, 1.1]),
    norm_ref_point=False,
    ideal=approx_idel,
    nadir=approx_nadir
)

# 根据算法求得的最优解计算Hypervolume指标
hv = [metric.do(_F) for _F in hist_F]

# 绘制收敛分析结果
plt.figure(figsize=(7, 5))
plt.plot(n_evals, hv,  color='black', lw=0.7, label="Avg. CV of Pop")
plt.scatter(n_evals, hv,  facecolor="none", edgecolor='black', marker="p")
plt.title("Convergence")
plt.xlabel("Function Evaluations")
plt.ylabel("Hypervolume")
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

代码执行结果如下图所示：

🚀 说明：
随着维数的增加，Hypervolume的计算成本会越来越高。精确的超体积可以有效地计算2个和3个目标。对于更高维度，目前Pymoo中还没有实现。

3.2 基于Running Metric的收敛性分析

当真实的Pareto前沿未知时，还可以使用近几年提出的性能评价指标Running Metric。Running Metric指标度量显示了不同代之间目标空间的差异，并使用算法的存活率来可视化改进。在Pymoo中，如果没有定义终止标准，也可以使用Running Metric作为多目标优化算法的终止条件。

比如，下面的分析表明，算法从第4代到第5代有了显著改进。

from pymoo.util.running_metric import RunningMetricAnimation

running = RunningMetricAnimation(
    delta_gen=5,
    n_plots=3,
    key_press=False,
    do_show=True
)

for algorithm in res.history[:15]:
    running.update(algorithm)
1
2
3
4
5
6
7
8
9
10
11

代码执行结果如下图所示：

绘制直到最终填充的图显示，该算法似乎收敛速度较慢，并且只做了轻微改进，此时终止算法迭代。

from pymoo.util.running_metric import RunningMetric

running = RunningMetricAnimation(delta_gen=10,
                        n_plots=4,
                        key_press=False,
                        do_show=True)

for algorithm in res.history:
    running.update(algorithm)
1
2
3
4
5
6
7
8
9

代码执行结果如下图所示：

📚参考文献

---

1. 🔍 (66条消息) Pymoo：使用多目标优化搜索解集的实现方法_嵌入式技术的博客-CSDN博客 ↩︎

2. 🔍 (66条消息) Python：带约束的双目标优化问题的实现方法_嵌入式技术的博客-CSDN博客 ↩︎

3. 🔍 (66条消息) Pymoo：优化算法的性能指标(Performance Indicators)_嵌入式技术的博客-CSDN博客 ↩︎ ↩︎