Overview#

This notebook demonstrates almost everything you can do with Paref, from applying an existing MOO algorithm to a bbf to constructing your own MOO algorithm tailored to the problem, using a practical example.

Setup#

Let us imagine the following situation: We have a blackbox function (bbf) that we want to optimize multicriterially or, synonymously, multiobjectively (multi objective optimization, MOO). In other words, we are searching for the bbf Pareto front. The bbf can represent a machine for which we can set certain parameters (vector \(x \in D\), where D represents the domain of the problem, which we also call design space) and which then reacts measurably via certain target variables (vector \(y \in T\), where T represents the co-domain of the problem, which we also call target space). Or the bbf represents a simulation model, where we can also define specifications x for its design and obtain responses y through simulation.

In our example, we choose as bbf the mathematical test function Zitzler-Deb-Thiele N.1. In the further course of our example, however, we assume that we do not know this function (which makes it a bbf). We can only “test” the bbf qua function calls in places we choose.

This approach offers the following advantage: ZDT-1 is not only defined analytically, but we also know (analytically) its Pareto front. Thus, we can immediately compare our results with the expected results and see how well our new approach works.

What Paref offers#

In many cases, the general MOO problem is formulated as “find the Pareto front” and not specified more clearly. Many MOO algorithms are accordingly designed to do just that: They identify Pareto-optimal points (i.e. points on the Pareto front), but they do not guarantee the user any further properties (e.g. to identify co-domain corner points of the Pareto front). The package paref now allows exactly that: besides the general property “identified points are Pareto-optimal”, paref users can specify further properties of the Pareto front to be identified and construct MOO algorithms based on those defined properties. On a top level, this works as follows:

  • properties of Pareto points are reflected by Pareto reflections or, more generally, by sequences of such

  • Paref provides a library of such (partly customizable) Pareto reflections including a description of what properties are targeted

  • applying an MOO to a blackbox function and a sequence results in an algorithmic search for Pareto points satisfying those properties (with mathematical guarantees).

This example#

In this example we exploit almost all the functionality Paref provides, from implementing a blackbox function and apply an existing Paref MOO algorithm reflecting some properties all the way to construct our own problem specific MOO algorithm. Concrete we will do the following:

  1. Implement a blackbox function (ZDT-2)

  2. Analyze the optimization process

  3. Apply one of Parefs’ MOO algorithms to the blackbox function

  4. Analyze the results

  5. Apply an MOO algorithm to some Pareto reflection in order to target a certain property of Pareto points

  6. Construct your own problem tailored MOO by using a customized Pareto reflection

Define and implement a blackbox function#

We first need to define how many dimensions our domain (Design Space) and co-domain (Target Space) have. In our example (ZDT-1) we want to map 3 Design Space dimensions to 2 Target Space dimensions. Second, we must define how it can retrieve function values at selected points, i.e. we must declare calls to our bbf. Those four information, i.e.

  • assignment of vector of design variables to vector of target values (here given by ZDT-1) implemented in the __call__ method

  • dimension of design space (here 3)

  • dimension of target space (here 2)

  • design space definition (here \([0,1]^{3}\))

must be implemented in Parefs` BlackboxFunction interface.

In addition, a proper scaling of the objectives is necessary in order for optimization to work properly as explained in trouble-shooting. In our case, we do not need to scale the objectives as they are already approximately in the same range.

[ ]:
import numpy as np
from paref.blackbox_functions.design_space.bounds import Bounds
from paref.interfaces.moo_algorithms.blackbox_function import BlackboxFunction


class ZDT1(BlackboxFunction):
    def __call__(self, x: np.ndarray) -> np.ndarray:
        n = len(x)
        f1 = x[0]
        g = 1 + 9 / (n - 1) * np.sum(x[1:])
        h = 1 - (f1 / g) ** (1 / 2)
        f2 = g * h
        return np.array([f1, f2])

    @property
    def dimension_design_space(self) -> int:
        return 3

    @property
    def dimension_target_space(self) -> int:
        return 2

    @property
    def design_space(self) -> Bounds:
        return Bounds(upper_bounds=np.ones(self.dimension_design_space),
                      lower_bounds=np.zeros(self.dimension_design_space))

We then initialize an instance of this bbf.

[ ]:
bbf = ZDT1()

The Pareto front of ZDT-1 looks as follows:

[ ]:
# This code is exclusively to visualize our results

import plotly.graph_objects as go

pareto_points_of_bbf = [
    i * np.eye(1, bbf.dimension_design_space, 0)[0]
    for i in np.arange(0, 1, 0.01)
]
pareto_front_of_bbf = np.array([bbf(point) for point in pareto_points_of_bbf])
bbf.clear_evaluations()

data = [
    go.Scatter(x=pareto_front_of_bbf.T[0],
               y=pareto_front_of_bbf.T[1],
               name='Real Pareto front',
               line=dict(width=4)),
]
fig = go.Figure(data=data)
fig.update_layout(title="Pareto front of ZDT-1",
                  width=600,
                  height=600,
                  plot_bgcolor='rgba(0,0,0,0)',
                  legend=dict(
                      x=0.2,
                      y=0.9,
                  ))
fig.show()

Note that in general we are not able to calculate the Pareto front. In this particular case, we can and use this knowledge in order to validate our results.

The initial exploration of the target space#

Every MOO algorithm needs to explore the target space first. This is typically done by a Latin Hypercube Sampling (LHC) which is implemented within the BlackboxFunction interface.

As a rule of thumb, you should start with 30 LHC evaluations.

[ ]:
bbf.perform_lhc(30)

Analyze the model fitness#

Paref strongly relies on the capability of the underlying model (GPR) to accurately approximate the bbf. Unfortunately, there exists no metric to measure the quality of the model in general but rather approximations. However, we can at least ensure that the training process of the model is successful and analyze the intrinsic uncertainty of the model.

[ ]:
from paref.express.info import Info

info = Info(bbf, training_iter=2000)  # 2000 is the default number of training iterations
[ ]:
info.model_fitness

We see that the loss of the GPRs converged and looks convex. This is how it should be. See trouble-shooting for more details on how to ensure the quality of the model and the resulting optimization.

Apply a MOO algorithm to blackbox function#

Applying an MOO algorithm in Paref is simply given by calling the algorithm to

  • a blackbox function (implemented in the BlackboxFunction interface) and

  • a stopping criterion indicating when the algorithm should terminate.

There are essentially two types of MOO algorithms implemented in Paref:

  • generic MOO algorithms (mainly minimizers) which are not tailored to some user defined properties but used when constructing a tailored MOO algorithm from a (sequence of) Pareto reflection(s)

  • tailored MOO algorithms which target certain properties of Pareto points

Let’s make it concrete for our example: From our bbf we want to identify Pareto-optimal solutions that represent the co-domain corners of the Pareto front. This property gives us an answer to the question “in which target area are the Pareto-optimal trade-offs of our bbf?”. This is a quite typical questions that arise in the context of a bbf MOO. Let’s get down to work and see how all these considerations transfer into code….

First, with initialize some stopping criteria. In our case, we choose the MaxIterationsReached which tells the algorithm to stop after a defined number max_iterations of iterations is reached. For an initial MOO, it is wise to grant the minimum number of iterations first. After analyzing the results, we can increase the number of iterations if necessary. Here, we grant one evaluation per component (i.e. per corner).

[ ]:
from paref.moo_algorithms.stopping_criteria.max_iterations_reached import MaxIterationsReached

stopping_criteria = MaxIterationsReached(max_iterations=2)

Next, we initialize the MOO algorithm targeting our desired properties of Pareto points: the FindEdgePoints algorithm.

[ ]:
from paref.moo_algorithms.multi_dimensional.find_edge_points import FindEdgePoints

moo_find_edge_points = FindEdgePoints()

At last, we simply apply the algorithm to the blackbox function and the stopping criteria.

[ ]:
moo_find_edge_points(blackbox_function=bbf,
                     stopping_criteria=stopping_criteria)

We access the so found best fitting Pareto points by calling the best_fits attribute of the algorithm.

[ ]:
moo_find_edge_points.best_fits
[ ]:
# This code is exclusively to visualize our results
data = [
    go.Scatter(x=pareto_front_of_bbf.T[0],
               y=pareto_front_of_bbf.T[1],
               name='Real Pareto front',
               line=dict(width=4)),
    go.Scatter(x=moo_find_edge_points.best_fits.T[0],
               y=moo_find_edge_points.best_fits.T[1],
               name='Determined Pareto points',
               mode="markers",
               marker=dict(size=10)),
]
fig = go.Figure(data=data)
fig.update_layout(title="Pareto front of ZDT-1",
                  width=600,
                  height=600,
                  plot_bgcolor='rgba(0,0,0,0)',
                  legend=dict(
                      x=0.2,
                      y=0.9,
                  ))
fig.show()

Analyze the results#

After (or also before) any MOO, we should analyze the results in order to guide further steps. For this task, Paref provides the Info class.

[ ]:
%%capture
info.update()

A central question is if the objectives are in fact conflicting (i.e. if there exists no global optimum) or, more generally, what the dimensionality of the Pareto front is.

[ ]:
info.topology

After convincing ourselves that the objectives are conflicting (so it is worthwhile to perform a MOO) we can check what minimal values in each objective we can expect.

[ ]:
info.minima

Sometimes it is hard to figure out what Pareto points we should target. In that case, we can simply ask Paref to suggest Pareto points.

[ ]:
info.suggestion_pareto_points

Identify Pareto points with defined properties by using Pareto reflections#

If an additional property of Pareto points (reflected by some Pareto reflection) is targeted, we may simply call the apply_to_sequence method of any MOO algorithm to the blackbox function, stopping criteria and the corresponding Pareto reflection in order to include this property.

Let’s make it concrete for our example: From our bbf ZDT-1 we want to determine the edge Pareto points of the Pareto front but restricted to a certain area. Here, we choose the area as \((-\infty,0.5]\times (-\infty,10]\), i.e. we demand the maximal \(y_1\) resp. \(y_2\) value of Pareto points to be \(0.5\) resp. \(10\).

This is a quite typical constrained that arises in the context of a bbf MOO.

Paref provides an implementation of a Pareto reflection corresponding to that property: the RestrictByPoint Pareto reflection.

[ ]:
from paref.pareto_reflections.restrict_by_point import RestrictByPoint

restricting_point = np.array([0.5, 10])
restrict_by_point = RestrictByPoint(restricting_point=restricting_point,
                                    nadir=10 * np.ones(2))

Now, we simply apply the FindEdgePoint algorithm to that Pareto reflection:

[ ]:
moo_find_edge_points_2 = FindEdgePoints()
moo_find_edge_points_2.apply_to_sequence(
    blackbox_function=bbf,
    sequence_pareto_reflections=restrict_by_point,
    stopping_criteria=MaxIterationsReached(2))
[ ]:
# This code is exclusively to visualize our results
area = np.array([[restricting_point[0], 0], restricting_point])
data = [
    go.Scatter(x=pareto_front_of_bbf.T[0],
               y=pareto_front_of_bbf.T[1],
               name='Real Pareto front',
               line=dict(width=4)),
    go.Scatter(x=bbf.y[-2:].T[0],
               y=bbf.y[-2:].T[1],
               name='Determined Pareto points',
               mode="markers",
               marker=dict(size=10)),
    go.Scatter(x=area.T[0],
               y=area.T[1],
               fill='tozerox',
               mode='none',
               fillcolor='rgba(255, 0, 0, 0.4)',
               name='Allowed area'),
]
fig = go.Figure(data=data)
fig.update_layout(title="Pareto front of ZDT-1",
                  width=600,
                  height=600,
                  plot_bgcolor='rgba(0,0,0,0)',
                  legend=dict(
                      x=0.5,
                      y=0.9,
                  ))
fig.show()

Customize Pareto reflections to user defined properties#

Imagine we prioritize Pareto points where the values of both objectives are close to each other. In other words, we want Pareto points \(y\) which minimize \(|y_1 - y_2|\). Clearly not every point which minimizes \(|y_1 - y_2|\) is Pareto optimal. However, for the purpose of finding Pareto points which minimize some function \(g\), there exists a universal Pareto reflection: the MinGParetoReflection. We simply need to specify \(g\) and the Pareto reflection does the rest.

[ ]:
from typing import Callable
from paref.pareto_reflections.minimize_g import MinGParetoReflection


class MinimizeDifferenceOfObjectives(MinGParetoReflection):
    @property
    def g(self) -> Callable:
        return lambda y: np.abs(y[0] - y[1])

    @property
    def dimension_domain(self) -> int:
        return 2


minimize_difference_of_objectives = MinimizeDifferenceOfObjectives(bbf)

In order to apply some MOO to that reflection (i.e. to target the ‘’minimize distance of objectives’’ property), we choose some generic MOO algorithm (i.e. some MOO which is not tailored to some properties) which can handle the MinimizeDifferenceOfObjectives reflection. The codomain dimension of the Pareto reflection (dimension_codomain property) must be a supported target space dimensions of the MOO (supported_target_space_dimensions) which is 1 in our case.

[ ]:
minimize_difference_of_objectives.dimension_codomain

In this example, we choose the DifferentialEvolutionMinimizer. This MOO exploits the fact that the underlying bbf is very cheap to sample (for a generic MOO algorithm which is tailored to expensive bbf see for example the GPRMinimizer) and yields typically better results but needs much (i.e. thousands of) more evaluations of the blackbox function.

[ ]:
from paref.moo_algorithms.minimizer.differential_evolution_minimizer import DifferentialEvolutionMinimizer

generic_moo = DifferentialEvolutionMinimizer()

Notice that the codomain dimension of the Pareto reflection is supported by the MOO:

[ ]:
print(
    f"Codomain dimension is supported: {minimize_difference_of_objectives.dimension_codomain in generic_moo.supported_codomain_dimensions}"
)

Again, we simply apply the MOO to the Pareto reflection by calling its apply_to_sequence method:

[ ]:
bbf.clear_evaluations()
generic_moo.apply_to_sequence(blackbox_function=bbf,
                              stopping_criteria=MaxIterationsReached(2),
                              sequence_pareto_reflections=minimize_difference_of_objectives)

Note how the algorithm perfectly found the (only) Pareto point which minimizes the difference of the objectives, i.e. has equal values for both objectives:

[ ]:
generic_moo.best_fits
[ ]:
# This code is exclusively to visualize our results
data = [
    go.Scatter(x=pareto_front_of_bbf.T[0],
               y=pareto_front_of_bbf.T[1],
               name='Real Pareto front',
               line=dict(width=4)),
    go.Scatter(x=bbf.pareto_front.T[0],
               y=bbf.pareto_front.T[1],
               name='Pareto point minimizing the distance of objectives',
               mode="markers",
               marker=dict(size=10)),
]
fig = go.Figure(data=data)
fig.update_layout(title="Pareto front of ZDT-1",
                  width=600,
                  height=600,
                  plot_bgcolor='rgba(0,0,0,0)',
                  legend=dict(
                      x=0.2,
                      y=0.9,
                  ))
fig.show()