Complete tutorial

In this section one can learn how to solve single-, multi-, and many-objective optimization problems using Differential Evolution based algorithms [3]. For more details about the algorithms used, please refer to the Algorithms section.

Let us start by importing python utilities for numerical operations, plotting, and optimization.

# Useful packages
import numpy as np
import matplotlib.pyplot as plt

# Problem classes
from pymoo.core.problem import Problem, ElementwiseProblem
from pymoo.problems.functional import FunctionalProblem
from pymoo.optimize import minimize

# Algorithmss
from pymoode.algorithms import DE, GDE3, NSDE, NSDER
from pymoo.algorithms.moo.nsga3 import NSGA3

# Reference directions and test problems
from pymoo.problems import get_problem
from pymoo.util.ref_dirs import get_reference_directions

# Survival operators
from pymoode.survival import RankAndCrowding, ConstrRankAndCrowding

# Termination
from pymoo.termination.default import DefaultSingleObjectiveTermination, DefaultMultiObjectiveTermination

# Performance
from pymoode.performance import SpacingIndicator

Defining a problem

There are three ways of defining a problem in pymoo:

1. Functional

2. Elementwise

3. Population-wise

Rastrigin

In this example, let us start by implementing the Rastrigin test function, described by the equations below.

\[\begin{split}\begin{align} \text{min} \; \; & f(\boldsymbol{x}) = An + \sum_{i=1}^{n}[x_i^2 - A \, \mathrm{cos}(2 \pi x_i)]\\ \text{s.t.} \; \; & -5.12 \leq x_i \leq 5.12 & \forall i \in \{ 1, 2 \} \end{align}\end{split}\]

Functional

In the functional implementation, each individual is evaluated separately and the user must parse objective functions and constraints as either list (in case more than one output is expected) or callables (in case just one objective or constraint is expected).

# Defining the objective function
def rastrigin(x):
    return np.sum(x * x - 10 * np.cos(2 * np.pi * x)) + 10 * np.size(x)


# Functional
functional_rastrigin = FunctionalProblem(
    n_var=2, objs=rastrigin,
    xl=np.full(2, -5.12), xu=np.full(2, 5.12),
    constr_ieq=[]
)

Elementwise

The Elementwise definition also evaluates individuals separately, but has some advantages of returning all objectives and constraints in the same method. Therefore it is possible to avoid unnecessary calculations performed in more than one objective or constraint.

class ElementwiseRasrigin(ElementwiseProblem):

    def __init__(self):

        xl = np.full(2, -5.12)
        xu = np.full(2, 5.12)

        super().__init__(n_var=2, n_obj=1, n_ieq_constr=0, xl=xl, xu=xu)

    def _evaluate(self, x, out, *args, **kwargs):
        out["F"] = rastrigin(x)

Population-wise

The Population-wise approach can be helpful if a vectorized numerical evaluation is possible, making solutions faster.

class Rastrigin(Problem):

    def __init__(self):

        xl = np.full(2, -5.12)
        xu = np.full(2, 5.12)

        super().__init__(n_var=2, n_obj=1, n_ieq_constr=0, xl=xl, xu=xu)

    def _evaluate(self, x, out, *args, **kwargs):

        out["F"] = np.sum(x * x - 10 * np.cos(2 * np.pi * x), axis=1) \
            + 10 * x.shape[1]

Single-objective DE

Let us use the Rastrigin problem as the example for single-objective DE [3].

de = DE(pop_size=30, variant="DE/rand/1/bin", F=(0.3, 1.0), CR=0.5)

NGEN = 100
SEED = 12

termination_single = DefaultSingleObjectiveTermination(
    xtol=1e-8,
    cvtol=0.0,
    ftol=1e-8,
    period=20,
    n_max_gen=NGEN,
)

rastrigin = Rastrigin()

res_de_problem_1 = minimize(
    rastrigin,
    de,
    termination=termination_single,
    seed=SEED,
    save_history=True,
    verbose=False
)

print(f"f: {(res_de_problem_1.F[0]):.4f}")
print(f"x_opt: {(res_de_problem_1.X)}")

f: 0.0000
x_opt: [ 1.90172541e-04 -2.75540458e-05]

# Creating a countouf plot
x1 = np.linspace(-5.12, 5.12, 100)
x2 = np.linspace(-5.12, 5.12, 100)
X1, X2 = np.meshgrid(x1, x2)
X = np.vstack((X1.flatten().reshape([1, -1]), X2.flatten().reshape([1, -1])))
F = np.array([rastrigin(x) for x in X.T])
F = F.reshape([100, 100])

fig, ax = plt.subplots(figsize=[6, 5], dpi=70)

fig.patch.set_facecolor('white')

surf = ax.contourf(X1, X2, F, levels=40)
fig.colorbar(surf, shrink=1.0, aspect=20)

ax.scatter(
    res_de_problem_1.history[9].pop.get("X")[:, 0],
    res_de_problem_1.history[9].pop.get("X")[:, 1],
    marker="x", color="yellow", label="10th gen",
)

ax.scatter(
    res_de_problem_1.history[-1].pop.get("X")[:, 0],
    res_de_problem_1.history[-1].pop.get("X")[:, 1],
    marker="x", color="red", label="Last gen",
)
ax.legend(framealpha=1.0)

ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")

plt.tight_layout()

Multi-objective problem

In this section, it will be introduced a problem with two objectives and two inequality constraints. You will notice that, in the object oriented approach, the objective functions are stored in the key "F" of the output dictionary and the inequalities are stored in the key "G".

In case there are equality constraints, one could store them in the key "H" of the output dictionary. However, I strongly advise to try and remove degrees of freedom of the problem by using equalities to reduce the number of variables instead of explicitly using equality constraints in DE.

Example functions

\[\begin{split}\begin{align} \text{min} \; \; & f_1(\boldsymbol{x}) = (x_1 - 0.5) ^ 2 + 0.7 x_1 x_2 + 1.2 (x_2 + 0.7) ^ 2\\ & f_2(\boldsymbol{x}) = (x_1 + 1.5) ^ 2 + 0.8 x_1 x_2 + 1.3 (x_2 - 1.7) ^ 2\\ \text{s.t.} \; \; & g_1(\boldsymbol{x}) \equiv x_1 ^ 2 + (x_2 - 1) ^ 2 - 9 \leq 0\\ & g_2(\boldsymbol{x}) \equiv - (x_1 + 0.5) ^ 2 - (x_2 - 1) ^ 2 + 2 \leq 0\\ & -5.0 \leq \boldsymbol{x} \leq 5.0 \end{align}\end{split}\]

class TwoObjExample(Problem):

    def __init__(self):

        xl = np.full(2, -5.0)
        xu = np.full(2, 5.0)

        super().__init__(n_var=2, n_obj=2, n_ieq_constr=2, xl=xl, xu=xu)

    def _evaluate(self, x, out, *args, **kwargs):

        F1 = (x[:, 0] - 0.5) ** 2 + 0.7 * x[:, 0] * \
            x[:, 1] + 1.2 * (x[:, 1] + 0.7) ** 2

        F2 = 1.1 * (x[:, 0] + 1.5) ** 2 + 0.8 * x[:, 0] * \
            x[:, 1] + 1.3 * (x[:, 1] - 1.7) ** 2

        out["F"] = np.column_stack([F1, F2])

        G1 = x[:, 0] ** 2 + (x[:, 1] - 1) ** 2 - 9
        G2 = - (x[:, 0] + 0.5) ** 2 - (x[:, 1] - 1) ** 2 + 2

        out["G"] = np.column_stack([G1, G2])

We will now instantiate two algorithms:

• GDE3 [6]

• NSDE

# Using a not so large value for CR due to the low dimensionality of the problem
gde3 = GDE3(pop_size=50, variant="DE/rand/1/bin", F=(0.0, 1.0), CR=0.7)
nsde = NSDE(pop_size=50, variant="DE/rand/1/bin", F=(0.0, 1.0), CR=0.7)

NGEN = 200

termination_multi = DefaultMultiObjectiveTermination(
    xtol=1e-8,
    cvtol=1e-8,
    ftol=1e-8,
    period=50,
    n_max_gen=NGEN,
)

res_gde3_problem_2 = minimize(
    TwoObjExample(),
    gde3,
    termination_multi,
    seed=SEED,
    save_history=True,
    verbose=False
)

fig, ax = plt.subplots(figsize=[6, 5], dpi=70)
ax.scatter(res_gde3_problem_2.F[:, 0],
           res_gde3_problem_2.F[:, 1], color="firebrick", label="GDE3")
ax.set_ylabel("$f_2$")
ax.set_xlabel("$f_1$")
ax.legend()
fig.tight_layout()
plt.show()

Pruning nondominated solutions

Although the results using the original GDE3 [6] already seemed great, they could have been improved by adopting the pruning strategy proposed by Kukkonen & Deb [11], which recursively re-calculates crowding distances as removes individuals from the population in the survival stage.

pruning_survival = RankAndCrowding(crowding_func="pcd")

gde3_ps = GDE3(pop_size=50, variant="DE/rand/1/bin",
               F=(0.0, 1.0), CR=0.7, survival=pruning_survival)

res_gde3_ps_problem_2 = minimize(
    TwoObjExample(),
    gde3_ps,
    ('n_gen', NGEN),
    seed=SEED,
    save_history=True,
    verbose=False
)

fig, ax = plt.subplots(1, 2, figsize=[12, 5], dpi=70)

ax[0].scatter(res_gde3_ps_problem_2.F[:, 0], res_gde3_ps_problem_2.F[:,
              1], color="indigo", label="GDE3-Pruning", marker="o")
ax[0].set_ylabel("$f_2$")
ax[0].set_xlabel("$f_1$")
ax[0].legend()

ax[1].scatter(res_gde3_problem_2.F[:, 0], res_gde3_problem_2.F[:, 1],
              color="firebrick", label="GDE3-Original", marker="o")
ax[1].set_ylabel("$f_2$")
ax[1].set_xlabel("$f_1$")
ax[1].legend()

fig.tight_layout()
plt.show()

Spacing indicator

By the figure above we could notice that the distribution of elements in the Pareto front has been improved by the pruning strategy. However, now we will use a quantitative metric to measure this improvement - the Spacing indicator.

It is calculated as it follows [12]:

\[\begin{align} S = \sqrt{ \frac{1}{|Q|} \sum_{i=1}^{Q} (d_{i} - \bar{d})^2 } \end{align}\]

In which, \(S\) is the spacing indicator, \(d_i\) is the distance of one element of index \(i\) to its nearest neighbor, \(\bar{d}\) is the mean of \(d_i\), and \(Q\) is the set of solutions.

By default the ‘cityblock’ distance metric is used, but any metric compatible with scipy.spatial.distance.pdist can be parsed when instantiating indicator via the keyword argument metric.

It might be interesting to normalize objectives, especially if they are in different scales. To do that, one should parse argument zero_to_one=True and provide either the pareto front pf or nadir and ideal points.

ideal = np.vstack(
    (res_gde3_problem_2.F, res_gde3_ps_problem_2.F)).min(axis=0) - 1e-3
nadir = np.vstack(
    (res_gde3_problem_2.F, res_gde3_ps_problem_2.F)).max(axis=0) + 1e-3

sp_p2 = SpacingIndicator(
    zero_to_one=True,
    ideal=ideal,
    nadir=nadir
)

print("SP GDE3", sp_p2.do(res_gde3_problem_2.F))
print("SP GDE3-PS", sp_p2.do(res_gde3_ps_problem_2.F))

SP GDE3 0.012024880243150753
SP GDE3-PS 0.004881671664926163

Many-objective problem

It is known that the original Rank and Crowding survival of NSGA-II [5] performs poorly in problems with more than two objectives [11] [7]. However, we can modify the algorithms to perform well in these problems by changing the survival operators.

For details about GDE3-MNN, refer to [7]; for details aboud NSDE-R, refer to [8]; and for detals about NSGA-III, refer to [9] and [10].

# Define the reference directions
ref_dirs = get_reference_directions("das-dennis", 3, n_partitions=15)

# Suggestion for NSGA-III
popsize = ref_dirs.shape[0] + ref_dirs.shape[0] % 4

gde3mnn = GDE3(pop_size=popsize, variant="DE/rand/1/bin", F=(0.0, 1.0), CR=0.2,
               survival=RankAndCrowding(crowding_func="mnn"))

nsder = NSDER(ref_dirs, pop_size=popsize,
              variant="DE/rand/1/bin", F=(0.0, 1.0), CR=0.5)

nsga3 = NSGA3(ref_dirs, pop_size=popsize)

DTLZ2

\[\begin{split}\begin{align} \text{min} \; \; & f_1(\boldsymbol{x}) = \mathrm{cos}(\pi / 2 x_1) \, \mathrm{cos}(\pi / 2 x_2) (1 + g(\boldsymbol{x}))\\ & f_2(\boldsymbol{x}) = \mathrm{cos}(\pi / 2 x_1) \, \mathrm{sin}(\pi / 2 x_2) (1 + g(\boldsymbol{x}))\\ & f_3(\boldsymbol{x}) = \mathrm{sin}(\pi / 2 x_1) (1 + g(\boldsymbol{x}))\\ \text{s.t.} \; \; & g(\boldsymbol{x}) = \sum_{i=1}^{3}(x_i - 0.5) ^ 2\\ & 0.0 \leq \boldsymbol{x} \leq 1.0 \end{align}\end{split}\]

dtlz2 = get_problem("dtlz2")

res_nsder_problem_3 = minimize(
    dtlz2,
    nsder,
    ('n_gen', 250),
    seed=SEED,
    save_history=True,
    verbose=False
)

res_gde3_problem_3 = minimize(
    dtlz2,
    gde3mnn,
    ('n_gen', 250),
    seed=SEED,
    save_history=True,
    verbose=False
)

res_nsga3_problem_3 = minimize(
    dtlz2,
    nsga3,
    ('n_gen', 250),
    seed=SEED,
    save_history=True,
    verbose=False
)

fig, ax = plt.subplots(figsize=[6, 5], dpi=70,
                       subplot_kw={'projection': '3d'})

ax.scatter(res_nsder_problem_3.F[:, 0], res_nsder_problem_3.F[:, 1], res_nsder_problem_3.F[:, 2],
           color="firebrick", label="NSDE-R", marker="o")

ax.scatter(res_gde3_problem_3.F[:, 0], res_gde3_problem_3.F[:, 1], res_gde3_problem_3.F[:, 2],
           color="navy", label="GDE3-MNN", marker="o")

ax.scatter(res_nsga3_problem_3.F[:, 0], res_nsga3_problem_3.F[:, 1], res_nsga3_problem_3.F[:, 2],
           color="grey", label="NSGA-III", marker="o")

ax.view_init(elev=30, azim=30)

ax.set_xlabel("$f_1$")
ax.set_ylabel("$f_2$")
ax.set_zlabel("$f_3$")
ax.legend()
fig.tight_layout()