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()
输出如下: