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$
# set values if testing
import os
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
# Ignore all warnings
import warnings
warnings.filterwarnings("ignore")
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
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.
from copy import deepcopy
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
generator.numerical_optimizer.n_restarts = NUM_RESTARTS
generator.n_monte_carlo_samples = N_MC_SAMPLES
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
Xopt ________________________________ Version: 0+untagged.1.ge872ea5 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 fixed_features: null gp_constructor: covar_modules: {} custom_noise_prior: null mean_modules: {} name: standard trainable_mean_keys: [] transform_inputs: true use_low_noise_prior: true log_transform_acquisition_function: false 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: null n_restarts: 20 name: LBFGS reference_point: s: 0.0 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: {} 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
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.calculate_hypervolume():.4}")
n_samples: 3 budget used: 0.01031 hypervolume: 0.0375
n_samples: 4 budget used: 0.01888 hypervolume: 0.0375
n_samples: 5 budget used: 0.02889 hypervolume: 0.2629
n_samples: 6 budget used: 0.03889 hypervolume: 0.5344
n_samples: 7 budget used: 0.06441 hypervolume: 0.7889
n_samples: 8 budget used: 0.1124 hypervolume: 0.9101
n_samples: 9 budget used: 0.2039 hypervolume: 1.136
n_samples: 10 budget used: 0.3941 hypervolume: 1.4
n_samples: 11 budget used: 0.8097 hypervolume: 1.751
n_samples: 12 budget used: 1.653 hypervolume: 2.143
n_samples: 13 budget used: 2.653 hypervolume: 2.25
n_samples: 14 budget used: 3.653 hypervolume: 2.25
n_samples: 15 budget used: 3.68 hypervolume: 2.25
n_samples: 16 budget used: 4.68 hypervolume: 2.25
n_samples: 17 budget used: 5.68 hypervolume: 2.25
n_samples: 18 budget used: 5.902 hypervolume: 2.25
n_samples: 19 budget used: 6.164 hypervolume: 2.25
n_samples: 20 budget used: 6.292 hypervolume: 2.25
n_samples: 21 budget used: 6.292 hypervolume: 2.25
n_samples: 22 budget used: 6.317 hypervolume: 2.25
n_samples: 23 budget used: 7.216 hypervolume: 2.25
n_samples: 24 budget used: 7.606 hypervolume: 2.25
n_samples: 25 budget used: 8.606 hypervolume: 2.25
n_samples: 26 budget used: 8.886 hypervolume: 2.25
n_samples: 27 budget used: 9.886 hypervolume: 2.25
n_samples: 28 budget used: 10.59 hypervolume: 2.25
Show resultsĀ¶
X.data
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.000059 | False |
1 | 0.750000 | 1.000000 | 0.100000 | dummy_constant | 0.750000 | 1.000000 | 0.626888 | 0.312500 | 0.000012 | False |
2 | 0.221312 | 1.671585 | 0.268254 | dummy_constant | 0.221312 | 1.671585 | 1.894185 | 1.450279 | 0.000032 | False |
3 | 1.731745 | 3.141590 | 0.256705 | dummy_constant | 1.731745 | 3.141590 | 11.889047 | 8.495194 | 0.000028 | False |
4 | 0.000000 | 0.846765 | 0.268271 | dummy_constant | 0.000000 | 0.846765 | -0.382990 | 0.370246 | 0.000027 | False |
5 | 0.000000 | 0.172064 | 0.268269 | dummy_constant | 0.000000 | 0.172064 | -1.070394 | 0.357542 | 0.000028 | False |
6 | 0.000000 | 0.000000 | 0.350619 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000025 | False |
7 | 0.333845 | 0.000000 | 0.419895 | dummy_constant | 0.333845 | 0.000000 | -0.988548 | 0.277607 | 0.000028 | False |
8 | 0.000000 | 0.000000 | 0.505057 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000026 | False |
9 | 0.000000 | 0.000000 | 0.622306 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000026 | False |
10 | 0.000000 | 0.000000 | 0.778154 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000024 | False |
11 | 0.000000 | 0.000000 | 0.952472 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000026 | False |
12 | 0.000000 | 0.000000 | 1.000000 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000026 | False |
13 | 0.000000 | 0.604637 | 1.000000 | dummy_constant | 0.000000 | 0.604637 | -0.734414 | 0.260949 | 0.000026 | False |
14 | 3.141590 | 0.000000 | 0.355001 | dummy_constant | 3.141590 | 0.000000 | 8.769588 | 7.227998 | 0.000027 | False |
15 | 0.595660 | 0.000000 | 1.000000 | dummy_constant | 0.595660 | 0.000000 | -0.745189 | 0.259151 | 0.000026 | False |
16 | 0.272524 | 0.000000 | 1.000000 | dummy_constant | 0.272524 | 0.000000 | -1.025731 | 0.301745 | 0.000026 | False |
17 | 0.911935 | 0.000000 | 0.650526 | dummy_constant | 0.911935 | 0.000000 | -0.268375 | 0.419690 | 0.000028 | False |
18 | 2.071884 | 1.084787 | 0.682534 | dummy_constant | 2.071884 | 1.084787 | 4.455836 | 2.812794 | 0.000028 | False |
19 | 0.000000 | 0.555851 | 0.555293 | dummy_constant | 0.000000 | 0.555851 | -0.791029 | 0.253119 | 0.000030 | False |
20 | 0.000000 | 0.000000 | 0.112890 | dummy_constant | 0.000000 | 0.000000 | -1.100000 | 0.500000 | 0.000027 | False |
21 | 1.728715 | 0.075189 | 0.346768 | dummy_constant | 1.728715 | 0.075189 | 1.917335 | 1.690206 | 0.000026 | False |
22 | 2.250353 | 0.714947 | 0.969946 | dummy_constant | 2.250353 | 0.714947 | 4.554438 | 3.109937 | 0.000027 | False |
23 | 2.368272 | 1.864756 | 0.764129 | dummy_constant | 2.368272 | 1.864756 | 8.117814 | 5.352998 | 0.000028 | False |
24 | 1.121965 | 0.000000 | 1.000000 | dummy_constant | 1.121965 | 0.000000 | 0.158806 | 0.636841 | 0.000027 | False |
25 | 0.000000 | 1.028563 | 0.695515 | dummy_constant | 0.000000 | 1.028563 | -0.042058 | 0.529379 | 0.000028 | False |
26 | 0.569697 | 0.480709 | 1.000000 | dummy_constant | 0.569697 | 0.480709 | -0.466045 | 0.005230 | 0.000029 | False |
27 | 1.399906 | 1.826007 | 0.903700 | dummy_constant | 1.399906 | 1.826007 | 4.344627 | 2.568124 | 0.000027 | 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).
import matplotlib.pyplot as plt
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Ā¶
ax = history.hist(["x1", "x2", "s"],bins=20)
history.plot(y=["x1", "x2", "s"])
<Axes: >
Plot the acqusisition functionĀ¶
Here we plot the acquisition function at a small set of fidelities $[0, 0.5, 1.0]$.
# 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.58072495 1.3505124 ]]
[<matplotlib.lines.Line2D at 0x7fe2b80b0430>]
# examine lengthscale of the first objective
list(model.models[0].named_parameters())
[('likelihood.noise_covar.raw_noise', Parameter containing: tensor([-25.3393], dtype=torch.float64, requires_grad=True)), ('mean_module.raw_constant', Parameter containing: tensor(-0.2050, dtype=torch.float64, requires_grad=True)), ('covar_module.raw_outputscale', Parameter containing: tensor(0.0106, dtype=torch.float64, requires_grad=True)), ('covar_module.base_kernel.raw_lengthscale', Parameter containing: tensor([[0.5430, 1.9576, 1.7260]], dtype=torch.float64, requires_grad=True))]