Pymoo学习 (5)：收敛性分析

1 引入

  收敛性分析可以帮助我们在再次解决同一个问题时，更好地决定终止标准。收敛性分析应当考虑两种情况：
  1）Pareto前沿未知；
  2）Pareto前沿通过解析推导而来，或者存在合理的近似值。
  一个双目标优化示意如https://inkiyinji.blog.csdn.net//article/details/125914306，其Pareto前沿和解如下：

  代码如下：

from pymoo.util.misc import stack

class MyTestProblem(MyProblem):

    def _calc_pareto_front(self, flatten=True, *args, **kwargs):
        f2 = lambda f1: ((f1/100) ** 0.5 - 1)**2
        F1_a, F1_b = np.linspace(1, 16, 300), np.linspace(36, 81, 300)
        F2_a, F2_b = f2(F1_a), f2(F1_b)

        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)

    def _calc_pareto_set(self, *args, **kwargs):
        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)


import matplotlib.pyplot as plt


problem = MyTestProblem()
pf_a, pf_b = problem.pareto_front(use_cache=False, flatten=False)
pf = problem.pareto_front(use_cache=False, flatten=True)

plt.figure(figsize=(7, 5))
plt.scatter(F[:, 0], F[:, 1], s=30, facecolors='none', edgecolors='b', label="Solutions")
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()

  示意中开启了一个历史记录选项save_history，这正是收敛性分析所需要的。当开启save_history后，其将记录每次迭代的算法对象并存储于Result对象中。该方法更占内存 (特别是对于多次迭代)，但具有可以事后分析任何算法相关变量的优点。
  为了进行分析，算法的一些中间步骤需要被考虑。它们已经被以下实现：
  1）Callback类存储了算法每次迭代的必要信息；
  2）当使用最小化方法时开启save_history以记录。
  此外使用Hypervolume和IGD来衡量算法的性能：

from pymoo.optimize import minimize

res = minimize(problem,
               algorithm,
               ("n_gen", 40),
               seed=1,
               save_history=True,
               verbose=False)

X, F = res.opt.get("X", "F")

hist = res.history
print(len(hist))

  从history可以导出我们所需要的信息：

n_evals = []             # 函数评估次数
hist_F = []              # 每次迭代的目标值
hist_cv = []             # 每次迭代的约束违反
hist_cv_avg = []         # 所有种群的平均约束违反

for algo in hist:

    # 存储评估数
    n_evals.append(algo.evaluator.n_eval)
    # 从算法中检索最优
    opt = algo.opt

    # 存储约束违反
    hist_cv.append(opt.get("CV").min())
    hist_cv_avg.append(algo.pop.get("CV").mean())

    # 过滤出可行解，并添加至目标空间
    feas = np.where(opt.get("feasible"))[0]
    hist_F.append(opt.get("F")[feas])

2 约束满足

  首先，可以按如下查询每一代可行解的数量：

k = np.where(np.array(hist_cv) <= 0.0)[0].min()
print(f"经过{
      n_evals[k]}次评估，第{
      k}代至少有1个可行解")

  输出如下：

经过40次评估，第0代至少有1个可行解

  因为这个问题的复杂性低，所以很快就找到了可行的解决方案。然而，对应于自己的优化问题可能完全不同，因此需要先行分析：

  代码如下：

# 分析最不可行的最优解而非总体时，将’hist_cv_avg‘替换为’hist_cv‘
vals = hist_cv_avg
k = np.where(np.array(vals) <= 0.0)[0].min()

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")
plt.axvline(n_evals[k], color="red", label="All Feasible", linestyle="--")
plt.title("Convergence")
plt.xlabel("Function Evaluations")
plt.ylabel("Constraint Violation")
plt.legend()
plt.show()

3 Pareto前沿未知

  如果Pareto前沿未知，则无法知道算法是否已经收敛到真正的最优值。当然也不是没有任何进一步的信息。因为依然可以看到算法在优化过程中何时取得了大部分进展，并确定迭代次数是否增减。此外，以下指标用于将两种算法相互比较：1）Hypervolume和2）IGD。

3.1 Hypervolume (HV)

  Hypervolume是多目标优化中的一个很重要的指标，是Pareto兼容的，需要基于预定义参考点和提供的解决方案之间的volume。因此，首先定义参考点ref_point，其应当大于Pareto前沿的最大值：

approx_ideal = F.min(axis=0)
approx_nadir = F.max(axis=0)

from pymoo.indicators.hv import Hypervolume

metric = Hypervolume(ref_point= np.array([1.1, 1.1]),
                     norm_ref_point=False,
                     zero_to_one=True,
                     ideal=approx_ideal,
                     nadir=approx_nadir)

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

  输出如下：

  Hypervolume可以针对2–3个目标进行精确有效的计算，然而Hypervolume的时间成本随着维数的增加而增高。因此对于更高的维度，一些使用近似方法，这在Pymoo中尚不可用。

3.2 运行指标

  当真正的Pareto前沿未知时，另一种分析运行的方法是新近的运行指标。 运行指标显示了每一代目标空间的差异。如果没有定义默认的终止标准，这个指标也被用于Pymoo来确定终止条件。
  例如，该分析表明算法从第4代到第5代存在显着改进：

  代码如下：

from pymoo.util.running_metric import RunningMetric

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

for algorithm in res.history[:15]:
    running.notify(algorithm)

  设置更多的迭代次数时：

  代码如下：

from pymoo.util.running_metric import RunningMetric

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

for algorithm in res.history:
    running.notify(algorithm)

4 Pareto前沿已知或者近似

4.1 IGD/GD/IGD+/GD+

  对于一个问题，其Pareto前沿可以手动提供，也可以直接在问题定义中实现，以动态分析运行。这里展示了一个使用算法历史作为附加后处理步骤的示例。
  对于现实世界问题，必须使用近似值。可以通过运行算法几次并从所有解集中提取非支配解来获得近似值。如果只有一次运行，另一种方法是使用获得的非支配解集作为近似值。但是，结果仅表明算法在收敛到最终集合方面的进展程度。
  代码如下：

from pymoo.indicators.igd import IGD

metric = IGD(pf, zero_to_one=True)

igd = [metric.do(_F) for _F in hist_F]

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

  输出如下：

  代码如下：

from pymoo.indicators.igd_plus import IGDPlus

metric = IGDPlus(pf, zero_to_one=True)

igd = [metric.do(_F) for _F in hist_F]

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")
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】https://pymoo.org/getting_started/part_4.html