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

from xopt import Xopt, Evaluator
from xopt.generators.bayesian import MultiFidelityGenerator
from xopt.resources.test_functions.tnk import evaluate_TNK, tnk_vocs
from xopt.vocs import get_feasibility_data

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 from xopt import Xopt, Evaluator from xopt.generators.bayesian import MultiFidelityGenerator from xopt.resources.test_functions.tnk import evaluate_TNK, tnk_vocs from xopt.vocs import get_feasibility_data 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())

/home/runner/work/Xopt/Xopt/.venv/lib/python3.12/site-packages/torch/jit/_script.py:1491: FutureWarning: `torch.jit.script` is deprecated. Please switch to `torch.compile` or `torch.export`.
  warnings.warn(

{'variables': {'x1': {'dtype': None, 'default_value': None, 'domain': [0.0, 3.14159], 'type': 'ContinuousVariable'}, 'x2': {'dtype': None, 'default_value': None, 'domain': [0.0, 3.14159], 'type': 'ContinuousVariable'}}, 'objectives': {'y1': {'dtype': None, 'type': 'MinimizeObjective'}, 'y2': {'dtype': None, 'type': 'MinimizeObjective'}}, 'constraints': {'c1': {'dtype': None, 'value': 0.0, 'type': 'GreaterThanConstraint'}, 'c2': {'dtype': None, 'value': 0.5, 'type': 'LessThanConstraint'}}, 'constants': {'a': {'dtype': None, 'value': 'dummy_constant', 'type': '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} + 0.001\) as an example from a real life multi-fidelity simulation problem.

In [2]:

my_vocs = deepcopy(tnk_vocs)
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)

# 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) 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) # 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: 0.1.dev1+g2bcc306fc
Data size: 2
Config as YAML:
data_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
    train_config: null
    train_kwargs: null
    train_method: lbfgs
    train_model: true
    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:
    discrete_max_batch_size: 2048
    discrete_max_choices: 4096
    max_iter: 1000
    max_time: 5.0
    mixed_max_discrete_configurations: 512
    n_restarts: 20
    name: LBFGS
  reference_point:
    s: 0.0
    y1: 1.5
    y2: 1.5
  returns_id: false
  supports_batch_generation: true
  supports_constraints: true
  supports_contextual_variables: true
  supports_discrete_variables: true
  supports_multi_objective: true
  supports_no_objective: true
  turbo_controller: null
  use_cuda: false
  use_pf_as_initial_points: false
  vocs:
    constants:
      a:
        dtype: null
        type: Constant
        value: dummy_constant
    constraints:
      c1:
        dtype: null
        type: GreaterThanConstraint
        value: 0.0
      c2:
        dtype: null
        type: LessThanConstraint
        value: 0.5
    objectives:
      s:
        dtype: null
        type: MaximizeObjective
      y1:
        dtype: null
        type: MinimizeObjective
      y2:
        dtype: null
        type: MinimizeObjective
    observables: {}
    variables:
      s:
        default_value: null
        domain:
        - 0.0
        - 1.0
        dtype: null
        type: ContinuousVariable
      x1:
        default_value: null
        domain:
        - 0.0
        - 3.14159
        dtype: null
        type: ContinuousVariable
      x2:
        default_value: null
        domain:
        - 0.0
        - 3.14159
        dtype: null
        type: ContinuousVariable
serialize_inline: false
serialize_torch: false
stopping_condition: null
strict: true
xopt_dump_file: null

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.003318 hypervolume: 0.03941

n_samples: 4 budget used: 0.005099 hypervolume: 0.03941

n_samples: 5 budget used: 0.007303 hypervolume: 0.0894

n_samples: 6 budget used: 0.01022 hypervolume: 0.0894

n_samples: 7 budget used: 0.01198 hypervolume: 0.0894

n_samples: 8 budget used: 0.01947 hypervolume: 0.0894

n_samples: 9 budget used: 0.02726 hypervolume: 0.0894

n_samples: 10 budget used: 0.03644 hypervolume: 0.0894

n_samples: 11 budget used: 0.04527 hypervolume: 0.1834

n_samples: 12 budget used: 0.05939 hypervolume: 0.2777

n_samples: 13 budget used: 0.07897 hypervolume: 0.31

n_samples: 14 budget used: 0.1157 hypervolume: 0.31

n_samples: 15 budget used: 0.1648 hypervolume: 0.31

n_samples: 16 budget used: 0.228 hypervolume: 0.4003

n_samples: 17 budget used: 0.3223 hypervolume: 0.5005

n_samples: 18 budget used: 0.4956 hypervolume: 0.5779

n_samples: 19 budget used: 0.8147 hypervolume: 0.6628

n_samples: 20 budget used: 1.565 hypervolume: 0.6628

n_samples: 21 budget used: 1.615 hypervolume: 0.7181

n_samples: 22 budget used: 2.059 hypervolume: 0.8413

n_samples: 23 budget used: 3.06 hypervolume: 1.042

n_samples: 24 budget used: 3.063 hypervolume: 1.042

n_samples: 25 budget used: 4.064 hypervolume: 1.136

n_samples: 26 budget used: 5.065 hypervolume: 1.178

n_samples: 27 budget used: 5.098 hypervolume: 1.178

n_samples: 28 budget used: 5.866 hypervolume: 1.178

n_samples: 29 budget used: 5.988 hypervolume: 1.178

n_samples: 30 budget used: 6.177 hypervolume: 1.178

n_samples: 31 budget used: 6.408 hypervolume: 1.178

n_samples: 32 budget used: 7.409 hypervolume: 1.204

n_samples: 33 budget used: 7.423 hypervolume: 1.204

n_samples: 34 budget used: 8.424 hypervolume: 1.256

n_samples: 35 budget used: 8.495 hypervolume: 1.256

n_samples: 36 budget used: 9.496 hypervolume: 1.272

n_samples: 37 budget used: 10.5 hypervolume: 1.272

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.009127	False
1	0.750000	1.000000	0.100000	dummy_constant	0.750000	1.000000	0.626888	0.312500	0.003226	False
2	0.506397	1.164691	0.023418	dummy_constant	0.506397	1.164691	0.516808	0.441855	0.005561	False
3	0.386093	0.532099	0.129469	dummy_constant	0.386093	0.532099	-0.486327	0.014005	0.008313	False
4	0.779748	0.685051	0.146518	dummy_constant	0.779748	0.685051	0.026071	0.112503	0.005799	False
5	0.000000	0.086858	0.167241	dummy_constant	0.000000	0.086858	-1.092456	0.420686	0.003137	False
6	0.000000	3.141590	0.128677	dummy_constant	0.000000	3.141590	8.769588	7.227998	0.000140	False
7	0.866554	0.254906	0.237109	dummy_constant	0.866554	0.254906	-0.170655	0.194433	0.002769	False
8	1.186424	0.030931	0.240134	dummy_constant	1.186424	0.030931	0.317129	0.691204	0.000155	False
9	2.688397	0.411289	0.253335	dummy_constant	2.688397	0.411289	6.472303	4.796953	0.000148	False
10	1.051967	0.184143	0.250164	dummy_constant	1.051967	0.184143	0.233814	0.404434	0.000162	False
11	0.333444	0.975350	0.289893	dummy_constant	0.333444	0.975350	0.009523	0.253699	0.002261	False
12	0.199202	1.087987	0.320217	dummy_constant	0.199202	1.087987	0.320429	0.436208	0.000150	False
13	0.938836	0.366818	0.385886	dummy_constant	0.938836	0.366818	-0.078842	0.210315	0.000153	False
14	0.034375	0.984985	0.420401	dummy_constant	0.034375	0.984985	-0.113446	0.452018	0.000146	False
15	0.896237	0.645817	0.452241	dummy_constant	0.896237	0.645817	0.304738	0.178266	0.000148	False
16	0.098658	1.046080	0.507676	dummy_constant	0.098658	1.046080	0.097395	0.459279	0.000149	False
17	0.791297	0.710109	0.605102	dummy_constant	0.791297	0.710109	0.065492	0.129000	0.000146	False
18	0.482288	0.966955	0.720908	dummy_constant	0.482288	0.966955	0.123984	0.218360	0.000149	False
19	0.211291	0.823062	0.920774	dummy_constant	0.211291	0.823062	-0.214133	0.187722	0.000150	False
20	1.048526	0.085857	0.423069	dummy_constant	1.048526	0.085857	0.080724	0.472394	0.000147	False
21	1.068627	0.165355	0.792541	dummy_constant	1.068627	0.165355	0.246728	0.435324	0.000143	False
22	0.121108	1.036867	1.000000	dummy_constant	0.121108	1.036867	0.118318	0.431785	0.000143	False
23	0.813909	0.297080	0.160909	dummy_constant	0.813909	0.297080	-0.326830	0.139716	0.000144	False
24	1.003222	0.312886	1.000000	dummy_constant	1.003222	0.312886	0.091910	0.288244	0.000158	False
25	0.668950	0.792882	1.000000	dummy_constant	0.668950	0.792882	0.054569	0.114324	0.000146	False
26	0.362418	0.410444	0.375394	dummy_constant	0.362418	0.410444	-0.754809	0.026949	0.000140	False
27	2.586916	3.123385	0.926739	dummy_constant	2.586916	3.123385	15.440469	11.237366	0.000141	False
28	0.110705	1.807171	0.547454	dummy_constant	0.110705	1.807171	2.222331	1.860247	0.000144	False
29	1.267023	0.088758	0.620669	dummy_constant	1.267023	0.088758	0.569569	0.757445	0.000140	False
30	0.724704	0.076764	0.657024	dummy_constant	0.724704	0.076764	-0.457169	0.229621	0.000141	False
31	0.856004	0.547455	1.000000	dummy_constant	0.856004	0.547455	0.127345	0.128991	0.000142	False
32	1.427204	0.065283	0.286012	dummy_constant	1.427204	0.065283	0.966748	1.048686	0.000141	False
33	1.033836	0.058036	1.000000	dummy_constant	1.033836	0.058036	0.009808	0.480312	0.000147	False
34	0.903354	0.701811	0.468305	dummy_constant	0.903354	0.701811	0.350063	0.203422	0.000142	False
35	0.471723	0.906227	1.000000	dummy_constant	0.471723	0.906227	0.026365	0.165820	0.000145	False
36	0.900365	0.424198	1.000000	dummy_constant	0.900365	0.424198	-0.081778	0.166038	0.000142	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, get_feasibility_data(tnk_vocs, 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, get_feasibility_data(tnk_vocs, 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 acquisition function¶

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

In [8]:

fidelities = [0.0, 0.5, 1.0]
for fidelity in fidelities:
    X.generator.visualize_model(
        variable_names=["x1", "x2"],
        reference_point={"s": fidelity},
    )

fidelities = [0.0, 0.5, 1.0] for fidelity in fidelities: X.generator.visualize_model( variable_names=["x1", "x2"], reference_point={"s": fidelity}, )

In [9]:

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

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

Out[9]:

[('likelihood.noise_covar.raw_noise',
  Parameter containing:
  tensor([-66.1367], requires_grad=True)),
 ('mean_module.raw_constant',
  Parameter containing:
  tensor(0.6496, requires_grad=True)),
 ('covar_module.raw_lengthscale',
  Parameter containing:
  tensor([[ 0.4601, 25.6598, 37.3171]], requires_grad=True))]