Pymoo learning (5): convergence analysis

1 introduce

   Convergence analysis can help us solve the same problem again , Better determine the termination criteria . Convergence analysis Two situations should be considered ：
  1）Pareto The frontier is unknown ;
  2）Pareto The frontier is derived analytically , Or there is a reasonable approximation .
   One Bi objective optimization Schematic as https://inkiyinji.blog.csdn.net//article/details/125914306, Its Pareto Frontier reconciliation is as follows ：

   The code is as follows ：

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

   schematic A history option is enabled in save_history, This is exactly what convergence analysis requires . When open save_history after , It will record the algorithm object of each iteration and store it in Result In the object . This method takes up more memory ( Especially for multiple iterations ), But it has the advantage of analyzing any algorithm related variables afterwards .
   For analysis , Some intermediate steps of the algorithm need to be considered . They have been implemented by ：
  1）Callback Class stores the necessary information for each iteration of the algorithm ;
  2） Turn on when using the minimize method save_history To record .
   In addition, use Hypervolume and IGD To measure the performance of the algorithm ：

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

   from history We can export the information we need ：

n_evals = []             #  Function evaluation times 
hist_F = []              #  The target value of each iteration 
hist_cv = []             #  The constraint of each iteration violates 
hist_cv_avg = []         #  The average constraint of all species violates 

for algo in hist:

    #  Number of storage evaluations 
    n_evals.append(algo.evaluator.n_eval)
    #  Retrieve the best from the algorithm 
    opt = algo.opt

    #  Storage constraint violation 
    hist_cv.append(opt.get("CV").min())
    hist_cv_avg.append(algo.pop.get("CV").mean())

    #  Filter out feasible solutions , And add to the target space 
    feas = np.where(opt.get("feasible"))[0]
    hist_F.append(opt.get("F")[feas])

2 constraint satisfaction

   First , The number of feasible solutions of each generation can be queried as follows ：

k = np.where(np.array(hist_cv) <= 0.0)[0].min()
print(f" after {
      n_evals[k]} The assessment , The first {
      k} Generation at least 1 There is a feasible solution ")

   Output is as follows ：

 after 40 The assessment , The first 0 Generation at least 1 There is a feasible solution 

   Because the complexity of this problem is low , So we quickly found a feasible solution . However , The optimization problem corresponding to yourself may be completely different , Therefore, it needs to be analyzed first ：

   The code is as follows ：

#  When analyzing the most infeasible optimal solution rather than the population , take ’hist_cv_avg‘ Replace with ’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 The frontier is unknown

   If Pareto The frontier is unknown , It is impossible to know whether the algorithm has converged to the real optimal value . Of course, it's not without any further information . Because you can still see when the algorithm has made most progress in the optimization process , And determine whether the number of iterations increases or decreases . Besides , The following metrics are used to compare the two algorithms ：1）Hypervolume and 2）IGD.

3.1 Hypervolume (HV)

  Hypervolume It is a very important index in multi-objective optimization , yes Pareto compatible Of , Need to be based on Predefined reference points and Solutions provided Between volume. therefore , First define the reference point ref_point, It should be greater than Pareto The maximum value of the leading edge ：

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

   Output is as follows ：

  Hypervolume Can target 2–3 Accurate and effective calculation of targets , However Hypervolume The time cost of increases with the increase of dimension . So for higher dimensions , Some use approximate methods , This is in Pymoo Not yet available in .

3.2 Operation index

   When the real Pareto When the frontier is unknown , Another method of analyzing operations is recent Operation index . The operational indicators show the differences in the target space of each generation . If no default termination criteria are defined , This indicator is also used Pymoo To determine the termination conditions .
   for example , The analysis shows that the algorithm starts from 4 Generation to the third 5 There are significant improvements in generation ：

   The code is as follows ：

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)

   When setting more iterations ：

   The code is as follows ：

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 The frontier is known or approximate

4.1 IGD/GD/IGD+/GD+

   For a question , Its Pareto The leading edge can be provided manually , It can also be implemented directly in the problem definition , Run with dynamic analysis . Here is an example of using algorithm history as an additional post-processing step .
   For real world problems , Approximate values must be used . The approximation can be obtained by running the algorithm several times and extracting non dominated solutions from all solution sets . If only run once , Another method is to use the obtained set of non dominated solutions as an approximation . however , The results only show the progress of the algorithm in converging to the final set .
   The code is as follows ：

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

   Output is as follows ：

   The code is as follows ：

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

   Output is as follows ：

reference

【1】https://pymoo.org/getting_started/part_4.html