Multiobjective optimization with Pareto front mapping

One way to do multiobjective optimization is with Pareto optimization, which explores the set of Pareto-efficient points. A point is Pareto-efficient if there are no other valid points that are better at every objective: it shows the “trade-off” between several objectives.

[1]:
from blop.utils import prepare_re_env

%run -i $prepare_re_env.__file__ --db-type=temp
[2]:
from blop import DOF, Objective, Agent
from blop.utils import functions
from blop.dofs import BrownianMotion

import numpy as np


def digestion(df):
    for index, entry in df.iterrows():
        x1, x2 = entry.x1, entry.x2

        df.loc[index, "f1"] = (x1 - 2) ** 2 + (x2 - 1) + 2
        df.loc[index, "f2"] = 9 * x1 - (x2 - 1) + 2
        df.loc[index, "c1"] = x1**2 + x2**2
        df.loc[index, "c2"] = x1 - 3 * x2 + 10

    return df


dofs = [
    DOF(name="x1", search_domain=(-20, 20)),
    DOF(name="x2", search_domain=(-20, 20)),
]


objectives = [
    Objective(name="f1", target="min"),
    Objective(name="f2", target="min"),
    Objective(name="c1", target=(-np.inf, 225)),
    Objective(name="c2", target=(-np.inf, 0)),
]

agent = Agent(
    dofs=dofs,
    objectives=objectives,
    digestion=digestion,
    db=db,
)

(uid,) = RE(agent.learn("qr", n=64))


Transient Scan ID: 1     Time: 2024-09-26 23:34:10
Persistent Unique Scan ID: '6de1285f-dec1-4c03-8f29-d73e112d8f39'
New stream: 'primary'
+-----------+------------+------------+------------+
|   seq_num |       time |         x1 |         x2 |
+-----------+------------+------------+------------+
|         1 | 23:34:10.2 |     -1.038 |     -0.475 |
|         2 | 23:34:10.2 |     -0.608 |      0.570 |
|         3 | 23:34:10.2 |      4.819 |      1.053 |
|         4 | 23:34:10.2 |      4.311 |     -0.950 |
|         5 | 23:34:10.2 |      7.591 |     -4.854 |
|         6 | 23:34:10.2 |      5.077 |     -6.787 |
|         7 | 23:34:10.2 |      1.795 |     -7.879 |
|         8 | 23:34:10.2 |      9.613 |    -10.723 |
|         9 | 23:34:10.2 |     14.223 |     -8.912 |
|        10 | 23:34:10.2 |     17.503 |     -5.007 |
|        11 | 23:34:10.2 |     15.146 |     -3.197 |
|        12 | 23:34:10.2 |     11.864 |     -2.104 |
|        13 | 23:34:10.2 |     12.293 |      2.164 |
|        14 | 23:34:10.2 |     15.888 |      3.258 |
|        15 | 23:34:10.2 |     18.323 |      5.101 |
|        16 | 23:34:10.2 |     14.731 |      9.008 |
|        17 | 23:34:10.2 |     16.367 |     12.409 |
|        18 | 23:34:10.2 |     18.802 |     19.232 |
|        19 | 23:34:10.2 |     13.022 |     15.326 |
|        20 | 23:34:10.2 |     10.584 |     13.502 |
|        21 | 23:34:10.2 |      6.279 |     17.569 |
|        22 | 23:34:10.2 |      0.496 |     16.477 |
|        23 | 23:34:10.2 |     -4.706 |     16.940 |
|        24 | 23:34:10.2 |     -8.299 |     18.346 |
|        25 | 23:34:10.2 |    -15.844 |     19.695 |
|        26 | 23:34:10.2 |    -12.249 |     16.103 |
|        27 | 23:34:10.2 |     -5.932 |     10.170 |
|        28 | 23:34:10.2 |     -2.338 |     13.763 |
|        29 | 23:34:10.2 |      3.090 |     14.539 |
|        30 | 23:34:10.2 |      8.871 |     10.635 |
|        31 | 23:34:10.2 |      8.411 |      4.960 |
|        32 | 23:34:10.2 |      5.820 |      6.836 |
|        33 | 23:34:10.2 |      2.225 |      7.932 |
|        34 | 23:34:10.2 |     -2.978 |      8.727 |
|        35 | 23:34:10.2 |     -6.391 |      4.164 |
|        36 | 23:34:10.2 |     -8.758 |      7.320 |
|        37 | 23:34:10.2 |    -10.540 |      9.492 |
|        38 | 23:34:10.2 |    -14.775 |     13.037 |
|        39 | 23:34:10.2 |    -18.367 |     11.633 |
|        40 | 23:34:10.2 |    -16.323 |      5.896 |
|        41 | 23:34:10.2 |    -18.846 |      2.776 |
|        42 | 23:34:10.2 |    -19.667 |     -2.723 |
|        43 | 23:34:10.2 |    -17.065 |     -5.793 |
|        44 | 23:34:10.2 |    -13.574 |     -1.318 |
|        45 | 23:34:10.2 |    -13.065 |      1.367 |
|        46 | 23:34:10.2 |     -7.134 |     -4.068 |
|        47 | 23:34:10.2 |     -9.578 |     -7.260 |
|        48 | 23:34:10.3 |    -10.969 |     -9.385 |
|        49 | 23:34:10.3 |    -14.345 |    -13.105 |
+-----------+------------+------------+------------+
|   seq_num |       time |         x1 |         x2 |
+-----------+------------+------------+------------+
|        50 | 23:34:10.3 |    -17.625 |    -11.697 |
|        51 | 23:34:10.3 |    -15.024 |    -19.787 |
|        52 | 23:34:10.3 |    -11.742 |    -16.192 |
|        53 | 23:34:10.3 |     -7.557 |    -18.402 |
|        54 | 23:34:10.3 |     -4.277 |    -16.995 |
|        55 | 23:34:10.3 |     -5.111 |    -10.269 |
|        56 | 23:34:10.3 |     -3.485 |     -8.666 |
|        57 | 23:34:10.3 |     -1.829 |    -13.864 |
|        58 | 23:34:10.3 |      1.004 |    -16.541 |
|        59 | 23:34:10.3 |      3.520 |    -14.631 |
|        60 | 23:34:10.3 |      7.099 |    -17.636 |
|        61 | 23:34:10.3 |     11.092 |    -13.558 |
|        62 | 23:34:10.3 |     13.451 |    -15.425 |
|        63 | 23:34:10.3 |     17.187 |    -12.464 |
|        64 | 23:34:10.3 |     19.545 |    -19.333 |
+-----------+------------+------------+------------+
generator list_scan ['6de1285f'] (scan num: 1)



We can plot our fitness and constraint objectives to see their models:

[3]:
agent.plot_objectives()
../_images/tutorials_pareto-fronts_4_0.png

We can plot the Pareto front (the set of all Pareto-efficient points), which shows the trade-off between the two fitnesses. The points in blue comprise the Pareto front, while the points in red are either not Pareto efficient or are invalidated by one of the constraints.

[4]:
agent.plot_pareto_front()
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/blop/plotting.py:556: UserWarning: The use of `x.T` on tensors of dimension other than 2 to reverse their shape is deprecated and it will throw an error in a future release. Consider `x.mT` to transpose batches of matrices or `x.permute(*torch.arange(x.ndim - 1, -1, -1))` to reverse the dimensions of a tensor. (Triggered internally at ../aten/src/ATen/native/TensorShape.cpp:3697.)
  ax.scatter(*y[~constraint].T[[i, j]], c="r", marker="x", label="invalid")
../_images/tutorials_pareto-fronts_6_1.png

We can explore the Pareto front by choosing a random point on the Pareto front and computing the expected improvement in the hypervolume of all fitness objectives with respect to that point (called the “reference point”). All this is done automatically with the qnehvi acquisition function:

[5]:
# this is broken now but is fixed in the next PR
# RE(agent.learn("qnehvi", n=4))