The Platypus framework provides us with a Python library to solve and analyze multi-objective problems conveniently. In this notebook, we use Platypus to solve a multi-objective optimization problem – the 3 objective version of DTLZ2 (Deb et al., 2002a) – using the evolutionary algorithm NSGA-II (Deb et al., 2002b). We generate runtime visualizations for snapshots of the algorithm to help build an understanding of how the evolutionary algorithm works. Specifically, we are going to look at the population at each generation, their ranks (a key parameter NSGA-II uses to select offspring), the parallel coordinates plot of the current population, and runtime indicators such as hypervolume and generational distance.
Setting up the problem
I created a Google Colab for this post. Since this post is equivalent to the Google Colab file, you may choose to look at either one. This project is intended to be used in a Jupyter Notebook environment with Python.
First, we import the Python libraries we need for the visualization.
import numpy as np
import math
import pandas as pd
from pandas.plotting import parallel_coordinates
import random
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
from matplotlib import animation, rc, rcParams
rc('animation', html='jshtml')
from mpl_toolkits.mplot3d import Axes3D
!pip install platypus-opt
from platypus import (NSGAII, NSGAIII, DTLZ2, Hypervolume, EpsilonBoxArchive,
Solution, GenerationalDistance, InvertedGenerationalDistance,
Hypervolume, EpsilonIndicator, Spacing)
Our goal is to visualize the population at each generation of the evolutionary algorithm. To do so, we utilize an interface provided by the Platypus library: the callback function.
At each iteration of the algorithm, the callback function (if initialized) is called. We define our callback function to store the current number of function evaluations (algorithm.nfe) and all the data points in the current population (algorithm.result). Each population has the type Solution defined in Platypus, which not only contains the values of the variables but also attributes specific to the solver, such as ranks and crowding distances. We will access some of these attributes during visualization.
We also define the frequency of NFEs that we want to store the values. Saving a snapshot of the algorithm at every single iteration may be too expensive or unnecessary. Hence, we set a lapse for some number of iterations, in which the callback function does nothing unless the last NFE is more than the frequency amount apart from this NFE.
In the following example, we set up the problem DTLZ2 and use the callback function to solve it using NSGAII.
#define the frequency
frequency = 200
# define the problem definition
problem = DTLZ2(3)
# instantiate the optimization algorithm
algorithm = NSGAII(problem, divisions_outer=12)
# define callback function
solutions_list = []
hyp = []
nfe = []
last_calc = 0
def DTLZ2_callback(algorithm):
global last_calc
if algorithm.nfe == last_calc + frequency:
last_calc = algorithm.nfe
nfe.append(algorithm.nfe)
solutions_list.append(algorithm.result);
# optimize the problem using 10,000 function evaluations
algorithm.run(10000,DTLZ2_callback)
In order to calculate the metrics, we first need to define our reference set. For DTLZ2(3), our points lie on the unit sphere in the first octant. Let’s randomly select 1000 such points and plot them.
# generate the reference set for 3D DTLZ2
reference_set = EpsilonBoxArchive([0.02, 0.02, 0.02])
for _ in range(1000):
solution = Solution(problem)
solution.variables = [random.uniform(0,1) if i < problem.nobjs-1 else 0.5 for i in range(problem.nvars)]
solution.evaluate()
reference_set.add(solution)
fig_ref = plt.figure()
ax_ref = fig_ref.add_subplot(projection="3d")
ax_ref.scatter(
[s.objectives[0] for s in reference_set],
[s.objectives[1] for s in reference_set],
[s.objectives[2] for s in reference_set],
)
Given the reference set, we can now calculate the indicators for each population across iterations. The Platypus library provides us with all the functions we need to calculate the following indicators: generational distance, hypervolume, epsilon indicator, and spacing.
We initialize them in a dictionary and iterate over all saved populations to calculate the indicators.
# calculate the indicators
indicators = {"gd" : GenerationalDistance(reference_set),
"hyp" : Hypervolume(reference_set),
"ei" : EpsilonIndicator(reference_set),
"sp" : Spacing()}
indicator_results = {index : [] for index in indicators}
for indicator in tqdm(indicator_results):
for solution in tqdm(solutions_list):
indicator_results[indicator] += [indicators[indicator].calculate(solution)]
Setting up the visualization
At this point, we have the data we need to perform runtime visualizations. We will utilize the animation.FuncAnimation function in matplotlib to create interactive animations in Jupyter Notebook (or Google Colab). The idea behind creating such animations is to first initialize a static figure, and then define an update function to let the FuncAnimation know how to visualize new data for each iteration.
We define drawframe() that does the following: (1) clear the axis, so that previous data points are wiped out; (2) draw the new data points; (3) reset the limits of data axes so that the axes are consistent across frames; (4) update new data for indicator axes.
def drawframe(n):
# clear axes
ax.cla()
ax_parallel.cla()
# save results
result = solutions_list[n]
crowding_distances = [s.crowding_distance if s.crowding_distance != math.inf else 0 for s in result]
ranks = [s.rank for s in result]
result = solutions_list[n]
points = {
'X': [s.objectives[0] for s in result],
'Y': [s.objectives[1] for s in result],
'Z': [s.objectives[2] for s in result],
'rank': [s.rank for s in result],
'tag' : ['tag' for s in result]
}
df = pd.DataFrame(points)
# update new data points
ax.scatter(points['X'], points['Y'], points['Z'],
c = ranks,
alpha = 0.5,
linestyle="", marker="o",
cmap=cmap,
vmax=max_rank,
vmin=0)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_zlim(zlim)
ax.set_title('Solve DTLZ2, NFE = ' + str(nfe[n]))
# update the parallel coordinates plot
parallel_coordinates(df[['X','Y','Z','tag']], 'tag', ax=ax_parallel)
# update indicator plots
for indicator in indicator_axes:
indicator_ax = indicator_axes[indicator]
indicator_ax.plot(nfe[:n],indicator_results[indicator][:n], c = 'r')
indicator_ax.set_xlim(left = min(nfe), right=max(nfe))
indicator_ax.set_ylim(bottom = min(indicator_results[indicator]), top = max(indicator_results[indicator]))
With the drawframe() function, we create a static figure that initializes the axes we will feed into the drawframe() function. The initialization does the following: (1) set up the subplots of the figure; (2) calculate the maximum ranks of all points in all populations to determine the color mapping; (3) load the points from the first iteration; (4) initialize the scatter plots, parallel coordinates plot, and indicator plots.
fig = plt.figure(figsize=(20,10), dpi = 70)
fig.suptitle("Runtime Visualization", fontsize=20)
fig.subplots_adjust(wspace=0.3, hspace=0.3)
ax = fig.add_subplot(2, 3, 1, projection="3d")
ax_parallel = fig.add_subplot(2,3,4)
indicator_axes = {"gd" : fig.add_subplot(2, 3, 2),
"hyp" : fig.add_subplot(2, 3, 3),
"ei" : fig.add_subplot(2, 3, 5),
"sp" : fig.add_subplot(2, 3, 6)}
indicator_names = {"gd" : "Generational Distance",
"hyp" : "Hypervolume",
"ei" : "Epsilon Indicator",
"sp" : "Spacing"}
# load the ranks of all points
all_rank = [s.rank for result in solutions_list for s in result]
max_rank = max(all_rank)
# define the colormap
cmap = plt.get_cmap('Accent', max_rank)
# load the points from the first iteration
result = solutions_list[0]
points = {
'X': [s.objectives[0] for s in result],
'Y': [s.objectives[1] for s in result],
'Z': [s.objectives[2] for s in result],
'rank': [s.rank for s in result],
'tag' : ['tag' for s in result]
}
df = pd.DataFrame(points)
# create the scatter plot
graph = ax.scatter(
points['X'], points['Y'], points['Z'],
c = points['rank'],
alpha = 0.5,
cmap=cmap,
linestyle="", marker="o",
vmax=max_rank,
vmin=0
)
# create the parallel coordinates plot
parallel_coordinates(df[['X','Y','Z','tag']], 'tag', ax=ax_parallel)
plt.colorbar(graph, label='Rank', pad = 0.2)
# save the dimensions for later use
xlim = ax.get_xlim()
ylim = ax.get_ylim()
zlim = ax.get_zlim()
title = ax.set_title("DTLZ2 Static Figure")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
# initialize subplots for each indicator
for indicator in indicator_axes:
indicator_axes[indicator].plot(nfe[:0],indicator_results[indicator][:0])
indicator_axes[indicator].set_title(indicator_names[indicator] + " vs NFE")
indicator_axes[indicator].set_xlabel("NFE")
indicator_axes[indicator].set_ylabel(indicator_names[indicator])
Now we are ready to create an animation. We use the FuncAnimation() function, the initialized figure, and the drawframe() function to create the animation.
In the Jupyter Notebook or Google Colab, you will be able to play the animation using the play button at the bottom of the image. By default, the animation will play in a loop. You may select “once” so that the animation freezes in the last frame. Important: Before re-generating the animation, be sure to re-run the previous initialization so that the figure is reset. Otherwise, you may see overlapping points/lines.
ani = animation.FuncAnimation(fig, drawframe, frames=len(nfe), interval=20, blit=False)
ani
This visualization shows how the algorithm progresses as NFE grows. For the set of solutions, we clearly see how they converge to the reference set. Moreover, more and more points have lower ranks, indicating they are getting closer to the Pareto Front (points on the Pareto Front have rank = 0). The parallel coordinates plot shows how our solutions get narrowed down and the tradeoffs we could make. Finally, the four indicator plots track the performance of our algorithm as NFE increases. The trajectory of generational distance, hypervolume, and epsilon indicator suggests convergence.
In conclusion, the project highlights the potential of the Platypus library in Python in providing valuable insights into the progress of evolutionary algorithms, not just their final outcomes. Through the use of NSGA-II as an illustrative example, we have demonstrated the ability to monitor the ranks of points across generations. In Dave’s post, the runtime visualizations revealed the changing probabilities of variators across iterations. These findings emphasize the power of incorporating dynamic techniques to gain a comprehensive understanding of the runtime behavior of MOEA algorithms. I hope this project opens doors to further explore, visualize, and analyze the dynamics of evolutionary algorithms.
References
[1] Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002a). Scalable multi-objective optimization test problems. In Proceedings of the 2002 Congress on Evolutionary Computation. CEC’02 (Cat. No. 02TH8600) (Vol. 1, pp. 825-830). IEEE.
[2] Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002b). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on evolutionary computation, 6(2), 182-197.
[3] Gold, D. (2020, May 6). Beyond Hypervolume: Dynamic visualization of MOEA runtime. Water Programming: A Collaborative Research Blog. https://waterprogramming.wordpress.com/2020/05/06/beyond-hypervolume-dynamic-visualization-of-moea-runtime/
[4] Python – How to create 3D scatter animations – Stack Overflow https://stackoverflow.com/questions/41602588/how-to-create-3d-scatter-animations
[5] GitHub – Project-Platypus/Platypus: A Free and Open Source Python Library for Multiobjective Optimization https://github.com/Project-Platypus/Platypus
(1)
(2)
is a diagonal covariance matrix with ![[\Sigma_y]_{\alpha,\alpha}=\sigma^2_{\alpha, y}](https://s0.wp.com/latex.php?latex=%5B%5CSigma_y%5D_%7B%5Calpha%2C%5Calpha%7D%3D%5Csigma%5E2_%7B%5Calpha%2C+y%7D&bg=ffffff&fg=444444&s=0&c=20201002)
for each feature 
.
as:
(3)
(4)
(5)
(6)
(7)
Water Programming: A Collaborative Research Blog 