Multi-fidelity Multi-objective Bayesian Optimization¶

Here we attempt to solve for the constrained Pareto front of the TNK multi-objective optimization problem using Multi-Fidelity Multi-Objective Bayesian optimization. For simplicity we assume that the objective and constraint functions at lower fidelities is exactly equal to the functions at higher fidelities (this is obviously not a requirement, although for the best results lower fidelity calculations should correlate with higher fidelity ones). The algorithm should learn this relationship and use information gathered at lower fidelities to gather samples to improve the hypervolume of the Pareto front at the maximum fidelity.

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 import MultiFidelityGenerator
from xopt.resources.test_functions.tnk import evaluate_TNK, tnk_vocs

import matplotlib.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
BUDGET = 0.02 if SMOKE_TEST else 10

evaluator = Evaluator(function=evaluate_TNK)
print(tnk_vocs.dict())

# 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 import MultiFidelityGenerator from xopt.resources.test_functions.tnk import evaluate_TNK, tnk_vocs import matplotlib.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 BUDGET = 0.02 if SMOKE_TEST else 10 evaluator = Evaluator(function=evaluate_TNK) print(tnk_vocs.dict())

{'variables': {'x1': (0.0, 3.14159), 'x2': (0.0, 3.14159)}, 'constraints': {'c1': ('GREATER_THAN', 0.0), 'c2': ('LESS_THAN', 0.5)}, 'objectives': {'y1': 'MINIMIZE', 'y2': 'MINIMIZE'}, 'constants': {'a': 'dummy_constant'}, 'observables': []}

Set up the Multi-Fidelity Multi-objective optimization algorithm¶

Here we create the Multi-Fidelity generator object which can solve both single and multi-objective optimization problems depending on the number of objectives in VOCS. We specify a cost function as a function of fidelity parameter $s=[0,1]$ as $C(s) = s^{3.5}$ as an example from a real life multi-fidelity simulation problem.

In [2]:

my_vocs = deepcopy(tnk_vocs)
my_vocs.constraints = {}
generator = MultiFidelityGenerator(vocs=my_vocs, reference_point={"y1": 1.5, "y2": 1.5})

# set cost function according to approximate scaling of laser plasma accelerator
# problem, see https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.5.013063
generator.cost_function = lambda s: s**3.5 + 0.001
generator.numerical_optimizer.n_restarts = NUM_RESTARTS
generator.n_monte_carlo_samples = N_MC_SAMPLES
generator.gp_constructor.use_low_noise_prior = True

X = Xopt(generator=generator, evaluator=evaluator, vocs=my_vocs)

# evaluate at some explicit initial points
X.evaluate_data(pd.DataFrame({"x1": [1.0, 0.75], "x2": [0.75, 1.0], "s": [0.0, 0.1]}))

X

my_vocs = deepcopy(tnk_vocs) my_vocs.constraints = {} generator = MultiFidelityGenerator(vocs=my_vocs, reference_point={"y1": 1.5, "y2": 1.5}) # set cost function according to approximate scaling of laser plasma accelerator # problem, see https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.5.013063 generator.cost_function = lambda s: s**3.5 + 0.001 generator.numerical_optimizer.n_restarts = NUM_RESTARTS generator.n_monte_carlo_samples = N_MC_SAMPLES generator.gp_constructor.use_low_noise_prior = True X = Xopt(generator=generator, evaluator=evaluator, vocs=my_vocs) # evaluate at some explicit initial points X.evaluate_data(pd.DataFrame({"x1": [1.0, 0.75], "x2": [0.75, 1.0], "s": [0.0, 0.1]})) X

Out[2]:

            Xopt
________________________________
Version: 2.6.9.dev35+g35be778bc.d20251218
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: 1
  vectorized: false
generator:
  computation_time: null
  custom_objective: null
  fixed_features: null
  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
  max_travel_distances: null
  model: null
  n_candidates: 1
  n_interpolate_points: null
  n_monte_carlo_samples: 128
  name: multi_fidelity
  numerical_optimizer:
    max_iter: 2000
    max_time: 5.0
    n_restarts: 20
    name: LBFGS
  reference_point:
    s: 0.0
    y1: 1.5
    y2: 1.5
  supports_batch_generation: true
  supports_constraints: true
  supports_multi_objective: true
  turbo_controller: null
  use_cuda: false
  use_pf_as_initial_points: false
serialize_inline: false
serialize_torch: false
stopping_condition: null
strict: true
vocs:
  constants:
    a: dummy_constant
  constraints: {}
  objectives:
    s: MAXIMIZE
    y1: MINIMIZE
    y2: MINIMIZE
  observables: []
  variables:
    s:
    - 0
    - 1
    x1:
    - 0.0
    - 3.14159
    x2:
    - 0.0
    - 3.14159

Run optimization routine¶

Instead of ending the optimization routine after an explict number of samples we end optimization once a given optimization budget has been exceeded. WARNING: This will slightly exceed the given budget

In [3]:

budget = BUDGET
while X.generator.calculate_total_cost() < budget:
    X.step()
    print(
        f"n_samples: {len(X.data)} "
        f"budget used: {X.generator.calculate_total_cost():.4} "
        f"hypervolume: {X.generator.get_pareto_front_and_hypervolume()[-1]:.4}"
    )

budget = BUDGET while X.generator.calculate_total_cost() < budget: X.step() print( f"n_samples: {len(X.data)} " f"budget used: {X.generator.calculate_total_cost():.4} " f"hypervolume: {X.generator.get_pareto_front_and_hypervolume()[-1]:.4}" )

n_samples: 3 budget used: 0.003317 hypervolume: 0.03843

n_samples: 4 budget used: 0.004639 hypervolume: 0.05377

n_samples: 5 budget used: 0.006095 hypervolume: 0.05377

n_samples: 6 budget used: 0.007822 hypervolume: 0.2564

n_samples: 7 budget used: 0.01206 hypervolume: 0.4375

n_samples: 8 budget used: 0.02425 hypervolume: 0.6233

n_samples: 9 budget used: 0.05301 hypervolume: 0.808

n_samples: 10 budget used: 0.1204 hypervolume: 1.037

n_samples: 11 budget used: 0.2739 hypervolume: 1.315

n_samples: 12 budget used: 0.6176 hypervolume: 1.657

n_samples: 13 budget used: 1.362 hypervolume: 1.657

n_samples: 14 budget used: 2.255 hypervolume: 2.178

n_samples: 15 budget used: 2.409 hypervolume: 2.178

n_samples: 16 budget used: 3.41 hypervolume: 2.25

n_samples: 17 budget used: 3.987 hypervolume: 2.25

n_samples: 18 budget used: 4.156 hypervolume: 2.25

n_samples: 19 budget used: 4.172 hypervolume: 2.25

n_samples: 20 budget used: 4.195 hypervolume: 2.25

n_samples: 21 budget used: 4.667 hypervolume: 2.25

n_samples: 22 budget used: 4.789 hypervolume: 2.25

n_samples: 23 budget used: 4.983 hypervolume: 2.25

n_samples: 24 budget used: 4.984 hypervolume: 2.25

n_samples: 25 budget used: 5.186 hypervolume: 2.25

n_samples: 26 budget used: 5.431 hypervolume: 2.25

n_samples: 27 budget used: 5.432 hypervolume: 2.25

n_samples: 28 budget used: 5.836 hypervolume: 2.25

n_samples: 29 budget used: 5.846 hypervolume: 2.25

n_samples: 30 budget used: 5.857 hypervolume: 2.25

n_samples: 31 budget used: 5.884 hypervolume: 2.25

n_samples: 32 budget used: 6.012 hypervolume: 2.25

n_samples: 33 budget used: 6.186 hypervolume: 2.25

n_samples: 34 budget used: 6.301 hypervolume: 2.25

n_samples: 35 budget used: 7.219 hypervolume: 2.25

n_samples: 36 budget used: 7.26 hypervolume: 2.25

n_samples: 37 budget used: 7.289 hypervolume: 2.25

n_samples: 38 budget used: 7.293 hypervolume: 2.25

n_samples: 39 budget used: 7.511 hypervolume: 2.25

n_samples: 40 budget used: 7.527 hypervolume: 2.25

n_samples: 41 budget used: 7.638 hypervolume: 2.25

n_samples: 42 budget used: 7.817 hypervolume: 2.25

n_samples: 43 budget used: 7.975 hypervolume: 2.25

n_samples: 44 budget used: 8.142 hypervolume: 2.25

n_samples: 45 budget used: 8.144 hypervolume: 2.25

n_samples: 46 budget used: 9.038 hypervolume: 2.25

n_samples: 47 budget used: 9.716 hypervolume: 2.25

n_samples: 48 budget used: 10.18 hypervolume: 2.25

Show results¶

In [4]:

X.data

X.data

Out[4]:

	x1	x2	s	a	y1	y2	c1	c2	xopt_runtime	xopt_error
0	1.000000	0.750000	0.000000	dummy_constant	1.000000	0.750000	0.626888	0.312500	0.000266	False
1	0.750000	1.000000	0.100000	dummy_constant	0.750000	1.000000	0.626888	0.312500	0.000156	False
2	0.549607	1.199892	0.015442	dummy_constant	0.549607	1.199892	0.658673	0.492310	0.000150	False
3	1.144046	0.574217	0.100550	dummy_constant	1.144046	0.574217	0.598605	0.420304	0.000147	False
4	1.457075	2.410703	0.111005	dummy_constant	1.457075	2.410703	7.009325	4.566776	0.000150	False
5	0.152230	0.000000	0.126844	dummy_constant	0.152230	0.000000	-1.076826	0.370944	0.000148	False
6	0.000000	0.000000	0.194446	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000159	False
7	0.000000	0.000000	0.277004	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000145	False
8	0.000000	0.000000	0.359122	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000146	False
9	0.000000	0.000000	0.460812	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000145	False
10	0.000000	0.000000	0.584245	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000145	False
11	0.000000	0.000000	0.736430	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000146	False
12	3.141590	0.000000	0.918612	dummy_constant	3.141590	0.000000	8.769588	7.227998	0.000146	False
13	0.000000	0.000000	0.968054	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000146	False
14	0.527150	2.199038	0.584256	dummy_constant	0.527150	2.199038	4.194876	2.887467	0.000146	False
15	0.000000	0.000000	1.000000	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000145	False
16	0.000000	0.000000	0.854417	dummy_constant	0.000000	0.000000	-1.100000	0.500000	0.000146	False
17	2.035110	2.265834	0.600428	dummy_constant	2.035110	2.265834	8.210246	5.474734	0.000149	False
18	0.883328	2.544945	0.301289	dummy_constant	0.883328	2.544945	6.197865	4.328741	0.000147	False
19	0.182946	0.313909	0.336176	dummy_constant	0.182946	0.313909	-0.812454	0.135153	0.000145	False
20	1.066613	2.833808	0.806437	dummy_constant	1.066613	2.833808	8.081519	5.767711	0.000160	False
21	0.723361	0.048578	0.547147	dummy_constant	0.723361	0.048578	-0.522147	0.253672	0.000146	False
22	0.184522	0.241559	0.624795	dummy_constant	0.184522	0.241559	-0.854620	0.166318	0.000145	False
23	0.159321	1.226539	0.101948	dummy_constant	0.159321	1.226539	0.577368	0.643922	0.000138	False
24	0.995295	2.988243	0.631712	dummy_constant	0.995295	2.988243	8.878355	6.436671	0.000146	False
25	0.023841	0.808521	0.668549	dummy_constant	0.023841	0.808521	-0.434807	0.321912	0.000147	False
26	0.103977	3.036370	0.024387	dummy_constant	0.103977	3.036370	8.144979	6.590005	0.000147	False
27	0.753670	1.605331	0.771608	dummy_constant	0.753670	1.605331	2.071245	1.286106	0.000147	False
28	0.506690	0.385579	0.253569	dummy_constant	0.506690	0.385579	-0.539145	0.013137	0.000146	False
29	2.905351	1.517576	0.270165	dummy_constant	2.905351	1.517576	9.728934	6.821174	0.000149	False
30	2.483472	2.107125	0.351656	dummy_constant	2.483472	2.107125	9.581708	6.517013	0.000145	False
31	2.999748	2.273629	0.554646	dummy_constant	2.999748	2.273629	13.225863	9.394500	0.000146	False
32	0.266069	2.568326	0.605809	dummy_constant	0.266069	2.568326	5.675167	4.332695	0.000145	False
33	1.283728	1.112682	0.537872	dummy_constant	1.283728	1.112682	1.844267	0.989609	0.000157	False
34	0.718909	0.018101	0.975565	dummy_constant	0.718909	0.018101	-0.574840	0.280148	0.000150	False
35	0.938993	2.522569	0.398042	dummy_constant	0.938993	2.522569	6.161506	4.283500	0.000147	False
36	0.988558	0.159909	0.361807	dummy_constant	0.988558	0.159909	0.086701	0.354351	0.000109	False
37	0.353082	1.745983	0.188894	dummy_constant	0.353082	1.745983	2.272994	1.574059	0.000150	False
38	2.335090	1.016289	0.645785	dummy_constant	2.335090	1.016289	5.389518	3.634110	0.000152	False
39	1.652884	0.455987	0.305576	dummy_constant	1.652884	0.455987	1.979401	1.331079	0.000149	False
40	3.125238	2.746085	0.532365	dummy_constant	3.125238	2.746085	16.256770	11.936774	0.000150	False
41	0.153711	1.411984	0.610382	dummy_constant	0.153711	1.411984	1.033667	0.951630	0.000147	False
42	2.392980	1.251243	0.588552	dummy_constant	2.392980	1.251243	6.277471	4.147740	0.000151	False
43	0.332954	0.378858	0.599020	dummy_constant	0.332954	0.378858	-0.797055	0.042580	0.000145	False
44	1.899950	0.407075	0.145160	dummy_constant	1.899950	0.407075	2.872763	1.968496	0.000147	False
45	2.793640	1.801898	0.968220	dummy_constant	2.793640	1.801898	10.147924	6.955725	0.000151	False
46	1.950171	0.939538	0.894286	dummy_constant	1.950171	0.939538	3.623752	2.296191	0.000146	False
47	2.692957	0.134869	0.800298	dummy_constant	2.692957	0.134869	6.200583	4.942381	0.000148	False

Plot results¶

Here we plot the resulting observations in input space, colored by feasibility (neglecting the fact that these data points are at varying fidelities).

In [5]:

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")

Plot path through input space¶

In [6]:

ax = history.hist(["x1", "x2", "s"], bins=20)

ax = history.hist(["x1", "x2", "s"], bins=20)

In [7]:

history.plot(y=["x1", "x2", "s"])

history.plot(y=["x1", "x2", "s"])

Out[7]:

<Axes: >

Plot the acqusisition function¶

Here we plot the acquisition function at a small set of fidelities $[0, 0.5, 1.0]$.

In [8]:

# plot the acquisition function

bounds = X.generator.vocs.bounds
model = X.generator.model

# create mesh over non-fidelity parameters
n = 50
x = torch.linspace(*bounds.T[1], n)
y = torch.linspace(*bounds.T[2], n)
xx, yy = torch.meshgrid(x, y)

# plot function(s) at a single fidelity parameter
fidelities = [0.0, 0.5, 1.0]
for fidelity in fidelities:
    pts = torch.hstack([ele.reshape(-1, 1) for ele in (xx, yy)]).double()
    pts = torch.cat((torch.ones(pts.shape[0], 1) * fidelity, pts), dim=-1)

    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()

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

        c = ax.pcolor(xxn, yyn, acq.reshape(n, n), cmap="Blues")
        fig.colorbar(c)
        ax.set_title(f"Acquisition function - s: {fidelity}")

        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, "+")


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 over non-fidelity parameters n = 50 x = torch.linspace(*bounds.T[1], n) y = torch.linspace(*bounds.T[2], n) xx, yy = torch.meshgrid(x, y) # plot function(s) at a single fidelity parameter fidelities = [0.0, 0.5, 1.0] for fidelity in fidelities: pts = torch.hstack([ele.reshape(-1, 1) for ele in (xx, yy)]).double() pts = torch.cat((torch.ones(pts.shape[0], 1) * fidelity, pts), dim=-1) 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() xxn, yyn = xx.numpy(), yy.numpy() c = ax.pcolor(xxn, yyn, acq.reshape(n, n), cmap="Blues") fig.colorbar(c) ax.set_title(f"Acquisition function - s: {fidelity}") 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, "+") 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.44128711 2.15751276]]

Out[8]:

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

In [9]:

# examine lengthscale of the first objective
list(model.models[0].named_parameters())

# examine lengthscale of the first objective list(model.models[0].named_parameters())

Out[9]:

[('likelihood.noise_covar.raw_noise',
  Parameter containing:
  tensor([-160.2607], dtype=torch.float64, requires_grad=True)),
 ('mean_module.raw_constant',
  Parameter containing:
  tensor(0.0498, dtype=torch.float64, requires_grad=True)),
 ('covar_module.raw_lengthscale',
  Parameter containing:
  tensor([[ 0.7135, 50.5201, 51.3931]], dtype=torch.float64, requires_grad=True))]