Multi-objective Bayesian Optimization¶

TNK function $n=2$ variables: $x_i \in [0, \pi], i=1,2$

Objectives:

• $f_i(x) = x_i$

Constraints:

• $g_1(x) = -x_1^2 -x_2^2 + 1 + 0.1 \cos\left(16 \arctan \frac{x_1}{x_2}\right) \le 0$
• $g_2(x) = (x_1 - 1/2)^2 + (x_2-1/2)^2 \le 0.5$

In [1]:

# set values if testing
import os

from copy import deepcopy

import pandas as pd
import numpy as np
import torch

from xopt import Xopt, Evaluator
from xopt.generators.bayesian.mggpo import MGGPOGenerator
from xopt.resources.test_functions.tnk import evaluate_TNK, tnk_vocs

from xopt.generators.bayesian.objectives import feasibility

from matplotlib import pyplot as plt

# Ignore all warnings
import warnings

warnings.filterwarnings("ignore")

SMOKE_TEST = os.environ.get("SMOKE_TEST")
N_MC_SAMPLES = 1 if SMOKE_TEST else 128
NUM_RESTARTS = 1 if SMOKE_TEST else 20

evaluator = Evaluator(function=evaluate_TNK)
evaluator.max_workers = 10

# test check options
vocs = deepcopy(tnk_vocs)
gen = MGGPOGenerator(vocs=vocs, reference_point={"y1": 1.5, "y2": 1.5})
gen.n_monte_carlo_samples = N_MC_SAMPLES
gen.numerical_optimizer.n_restarts = NUM_RESTARTS
X = Xopt(evaluator=evaluator, generator=gen, vocs=vocs)
X.evaluate_data(pd.DataFrame({"x1": [1.0, 0.75], "x2": [0.75, 1.0]}))

X

# set values if testing import os from copy import deepcopy import pandas as pd import numpy as np import torch from xopt import Xopt, Evaluator from xopt.generators.bayesian.mggpo import MGGPOGenerator from xopt.resources.test_functions.tnk import evaluate_TNK, tnk_vocs from xopt.generators.bayesian.objectives import feasibility from matplotlib import pyplot as plt # Ignore all warnings import warnings warnings.filterwarnings("ignore") SMOKE_TEST = os.environ.get("SMOKE_TEST") N_MC_SAMPLES = 1 if SMOKE_TEST else 128 NUM_RESTARTS = 1 if SMOKE_TEST else 20 evaluator = Evaluator(function=evaluate_TNK) evaluator.max_workers = 10 # test check options vocs = deepcopy(tnk_vocs) gen = MGGPOGenerator(vocs=vocs, reference_point={"y1": 1.5, "y2": 1.5}) gen.n_monte_carlo_samples = N_MC_SAMPLES gen.numerical_optimizer.n_restarts = NUM_RESTARTS X = Xopt(evaluator=evaluator, generator=gen, vocs=vocs) X.evaluate_data(pd.DataFrame({"x1": [1.0, 0.75], "x2": [0.75, 1.0]})) X

Out[1]:

            Xopt
________________________________
Version: 2.4.6.dev5+ga295b108.d20250107
Data size: 2
Config as YAML:
dump_file: null
evaluator:
  function: xopt.resources.test_functions.tnk.evaluate_TNK
  function_kwargs:
    raise_probability: 0
    random_sleep: 0
    sleep: 0
  max_workers: 10
  vectorized: false
generator:
  computation_time: null
  custom_objective: null
  fixed_features: null
  ga_generator:
    crossover_probability: 0.9
    mutation_probability: 1.0
    output_path: null
    population: null
    population_file: null
    population_size: 64
    supports_multi_objective: true
  gp_constructor:
    covar_modules: {}
    custom_noise_prior: null
    mean_modules: {}
    name: standard
    trainable_mean_keys: []
    transform_inputs: true
    use_cached_hyperparameters: false
    use_low_noise_prior: true
  log_transform_acquisition_function: false
  max_travel_distances: null
  memory_length: null
  model: null
  n_candidates: 1
  n_interpolate_points: null
  n_monte_carlo_samples: 128
  name: mggpo
  numerical_optimizer:
    max_iter: 2000
    max_time: null
    n_restarts: 20
    name: LBFGS
  population_size: 64
  reference_point:
    y1: 1.5
    y2: 1.5
  supports_batch_generation: true
  supports_multi_objective: true
  turbo_controller: null
  use_cuda: false
max_evaluations: null
serialize_inline: false
serialize_torch: false
strict: true
vocs:
  constants:
    a: dummy_constant
  constraints:
    c1:
    - GREATER_THAN
    - 0.0
    c2:
    - LESS_THAN
    - 0.5
  objectives:
    y1: MINIMIZE
    y2: MINIMIZE
  observables: []
  variables:
    x1:
    - 0.0
    - 3.14159
    x2:
    - 0.0
    - 3.14159

In [2]:

for i in range(10):
    print(i)
    X.step()

for i in range(10): print(i) X.step()

0

1

2

3

4

5

6

7

8

9

In [3]:

X.generator.data

X.generator.data

Out[3]:

	x1	x2	a	y1	y2	c1	c2	xopt_runtime	xopt_error
0	1.000000	0.750000	dummy_constant	1.000000	0.750000	0.626888	0.312500	0.000161	False
1	0.750000	1.000000	dummy_constant	0.750000	1.000000	0.626888	0.312500	0.000133	False
2	0.299629	1.596399	dummy_constant	0.299629	1.596399	1.736772	1.242239	0.000151	False
3	0.310642	1.444926	dummy_constant	0.310642	1.444926	1.281284	0.928742	0.000131	False
4	0.316052	1.437940	dummy_constant	0.316052	1.437940	1.262481	0.913568	0.000128	False
...	...	...	...	...	...	...	...	...	...
97	0.094153	1.040030	dummy_constant	0.094153	1.040030	0.077934	0.456344	0.000101	False
98	0.912766	0.395104	dummy_constant	0.912766	0.395104	-0.107568	0.181379	0.000123	False
99	0.115304	1.019347	dummy_constant	0.115304	1.019347	0.075296	0.417712	0.000121	False
100	0.859603	0.506720	dummy_constant	0.859603	0.506720	0.057655	0.129359	0.000121	False
101	0.374796	0.944791	dummy_constant	0.374796	0.944791	-0.064016	0.213515	0.000121	False

102 rows × 9 columns

plot results¶

In [4]:

fig, ax = plt.subplots()

theta = np.linspace(0, np.pi / 2)
r = np.sqrt(1 + 0.1 * np.cos(16 * theta))
x_1 = r * np.sin(theta)
x_2_lower = r * np.cos(theta)
x_2_upper = (0.5 - (x_1 - 0.5) ** 2) ** 0.5 + 0.5

z = np.zeros_like(x_1)

# ax2.plot(x_1, x_2_lower,'r')
ax.fill_between(x_1, z, x_2_lower, fc="white")
circle = plt.Circle(
    (0.5, 0.5), 0.5**0.5, color="r", alpha=0.25, zorder=0, label="Valid Region"
)
ax.add_patch(circle)
history = pd.concat(
    [X.data, tnk_vocs.feasibility_data(X.data)], axis=1, ignore_index=False
)

ax.plot(*history[["x1", "x2"]][history["feasible"]].to_numpy().T, ".C1")
ax.plot(*history[["x1", "x2"]][~history["feasible"]].to_numpy().T, ".C2")

ax.set_xlim(0, 3.14)
ax.set_ylim(0, 3.14)
ax.set_xlabel("x1")
ax.set_ylabel("x2")
ax.set_aspect("equal")

fig, ax = plt.subplots() theta = np.linspace(0, np.pi / 2) r = np.sqrt(1 + 0.1 * np.cos(16 * theta)) x_1 = r * np.sin(theta) x_2_lower = r * np.cos(theta) x_2_upper = (0.5 - (x_1 - 0.5) ** 2) ** 0.5 + 0.5 z = np.zeros_like(x_1) # ax2.plot(x_1, x_2_lower,'r') ax.fill_between(x_1, z, x_2_lower, fc="white") circle = plt.Circle( (0.5, 0.5), 0.5**0.5, color="r", alpha=0.25, zorder=0, label="Valid Region" ) ax.add_patch(circle) history = pd.concat( [X.data, tnk_vocs.feasibility_data(X.data)], axis=1, ignore_index=False ) ax.plot(*history[["x1", "x2"]][history["feasible"]].to_numpy().T, ".C1") ax.plot(*history[["x1", "x2"]][~history["feasible"]].to_numpy().T, ".C2") ax.set_xlim(0, 3.14) ax.set_ylim(0, 3.14) ax.set_xlabel("x1") ax.set_ylabel("x2") ax.set_aspect("equal")

In [5]:

# plot model predictions

data = X.data

bounds = X.generator.vocs.bounds
model = X.generator.train_model(X.generator.data)

# create mesh
n = 50
x = torch.linspace(*bounds.T[0], n)
y = torch.linspace(*bounds.T[1], n)
xx, yy = torch.meshgrid(x, y)
pts = torch.hstack([ele.reshape(-1, 1) for ele in (xx, yy)]).double()

xx, yy = xx.numpy(), yy.numpy()

outputs = X.generator.vocs.output_names
with torch.no_grad():
    post = model.posterior(pts)

    for i in range(len(vocs.output_names)):
        mean = post.mean[..., i]
        fig, ax = plt.subplots()
        ax.plot(*data[["x1", "x2"]].to_numpy().T, "+C1")
        c = ax.pcolor(
            xx, yy, mean.squeeze().reshape(n, n), cmap="seismic", vmin=-10.0, vmax=10.0
        )
        fig.colorbar(c)
        ax.set_title(f"Posterior mean: {outputs[i]}")

# plot model predictions data = X.data bounds = X.generator.vocs.bounds model = X.generator.train_model(X.generator.data) # create mesh n = 50 x = torch.linspace(*bounds.T[0], n) y = torch.linspace(*bounds.T[1], n) xx, yy = torch.meshgrid(x, y) pts = torch.hstack([ele.reshape(-1, 1) for ele in (xx, yy)]).double() xx, yy = xx.numpy(), yy.numpy() outputs = X.generator.vocs.output_names with torch.no_grad(): post = model.posterior(pts) for i in range(len(vocs.output_names)): mean = post.mean[..., i] fig, ax = plt.subplots() ax.plot(*data[["x1", "x2"]].to_numpy().T, "+C1") c = ax.pcolor( xx, yy, mean.squeeze().reshape(n, n), cmap="seismic", vmin=-10.0, vmax=10.0 ) fig.colorbar(c) ax.set_title(f"Posterior mean: {outputs[i]}")

In [6]:

# plot the acquisition function
bounds = X.generator.vocs.bounds
model = X.generator.model

# create mesh
n = 25
x = torch.linspace(*bounds.T[0], n)
y = torch.linspace(*bounds.T[1], n)
xx, yy = torch.meshgrid(x, y)
pts = torch.hstack([ele.reshape(-1, 1) for ele in (xx, yy)]).double()

xx, yy = xx.numpy(), yy.numpy()

acq_func = X.generator.get_acquisition(model)
with torch.no_grad():
    acq_pts = pts.unsqueeze(1)
    acq = acq_func(acq_pts)

    fig, ax = plt.subplots()
    c = ax.pcolor(xx, yy, acq.reshape(n, n), cmap="Blues")
    fig.colorbar(c)
    ax.set_title("Acquisition function")

    ax.plot(*history[["x1", "x2"]][history["feasible"]].to_numpy().T, ".C1")
    ax.plot(*history[["x1", "x2"]][~history["feasible"]].to_numpy().T, ".C2")

    ax.plot(*history[["x1", "x2"]].to_numpy()[-1].T, "+")

    feas = feasibility(pts.unsqueeze(1), model, tnk_vocs).flatten()

    fig2, ax2 = plt.subplots()
    c = ax2.pcolor(xx, yy, feas.reshape(n, n))
    fig2.colorbar(c)
    ax2.set_title("Feasible Region")

candidate = pd.DataFrame(X.generator.generate(1), index=[0])
print(candidate[["x1", "x2"]].to_numpy())
ax.plot(*candidate[["x1", "x2"]].to_numpy()[0], "o")

# plot the acquisition function bounds = X.generator.vocs.bounds model = X.generator.model # create mesh n = 25 x = torch.linspace(*bounds.T[0], n) y = torch.linspace(*bounds.T[1], n) xx, yy = torch.meshgrid(x, y) pts = torch.hstack([ele.reshape(-1, 1) for ele in (xx, yy)]).double() xx, yy = xx.numpy(), yy.numpy() acq_func = X.generator.get_acquisition(model) with torch.no_grad(): acq_pts = pts.unsqueeze(1) acq = acq_func(acq_pts) fig, ax = plt.subplots() c = ax.pcolor(xx, yy, acq.reshape(n, n), cmap="Blues") fig.colorbar(c) ax.set_title("Acquisition function") ax.plot(*history[["x1", "x2"]][history["feasible"]].to_numpy().T, ".C1") ax.plot(*history[["x1", "x2"]][~history["feasible"]].to_numpy().T, ".C2") ax.plot(*history[["x1", "x2"]].to_numpy()[-1].T, "+") feas = feasibility(pts.unsqueeze(1), model, tnk_vocs).flatten() fig2, ax2 = plt.subplots() c = ax2.pcolor(xx, yy, feas.reshape(n, n)) fig2.colorbar(c) ax2.set_title("Feasible Region") candidate = pd.DataFrame(X.generator.generate(1), index=[0]) print(candidate[["x1", "x2"]].to_numpy()) ax.plot(*candidate[["x1", "x2"]].to_numpy()[0], "o")

[[0.5945952  0.10581534]]

Out[6]:

[<matplotlib.lines.Line2D at 0x7f1be469a490>]