Constrained single-point Bayesian optimisation#

In this example, NUBO is used to maximise a function where the input space is bounded and constrained. The whole process is not too different from the unconstrained case. We only need to choose a different optimiser that allows the use of constraints when maximising the acquisition function UpperConfidenceBound. NUBO uses the SLSQP optimiser that can be provided with a dictionary or a tuple of dictionaries that specify one or multiple constraints. We specify two constraints to showcase the two different options: equality constraints and inequality constraints. Equality constraints require the constraint to be 0 while the result is non-negative for inequality constraints. Our first constraint {'type': 'ineq', 'fun': lambda x: 0.5 - x[0] - x[1]} is an inequality constraint and requires the sum of the first two inputs to be smaller or equal to 0.5. The second constraint {'type': 'eq', 'fun': lambda x: 1.2442 - x[3] - x[4] - x[5]} is an equality constraint specifying that the sum of the last three inputs needs to be equal to 1.2442. These constraints are very simple and in practice, much more complex constraints might be specified. The Hartmann6D synthetic test function acts as a substitute for a black-box objective function, such as an experiment or a simulation. The optimisation loop is run for 40 iterations and finds a solution close to the true optimum of -3.3224. Important: Generating initial input points with a Latin hypercube might not work for real problems as they will not consider the constraints but only the bounds. In these situations, other methods or selecting initial points by hand might be preferable. The purpose of this example is solely the demonstration of how NUBO handles constraints and constrained optimisation.

[1]:

import torch
from nubo.acquisition import ExpectedImprovement, UpperConfidenceBound
from nubo.models import GaussianProcess, fit_gp
from nubo.optimisation import single
from nubo.test_functions import Hartmann6D
from nubo.utils import gen_inputs
from gpytorch.likelihoods import GaussianLikelihood


# test function
func = Hartmann6D(minimise=False)
dims = 6

# specify bounds and constraints
bounds = torch.tensor([[0., 0., 0., 0., 0., 0.], [1., 1., 1., 1., 1., 1.]])
cons = ({'type': 'ineq', 'fun': lambda x: 0.5 - x[0] - x[1]},
        {'type': 'eq', 'fun': lambda x: 1.2442 - x[3] - x[4] - x[5]})

# training data
x_train = gen_inputs(num_points=dims*5,
                     num_dims=dims,
                     bounds=bounds)
y_train = func(x_train)

# Bayesian optimisation loop
iters = 40

for iter in range(iters):

    # specify Gaussian process
    likelihood = GaussianLikelihood()
    gp = GaussianProcess(x_train, y_train, likelihood=likelihood)

    # fit Gaussian process
    fit_gp(x_train, y_train, gp=gp, likelihood=likelihood, lr=0.1, steps=200)

    # specify acquisition function
    # acq = ExpectedImprovement(gp=gp, y_best=torch.max(y_train))
    acq = UpperConfidenceBound(gp=gp, beta=1.96**2)

    # optimise acquisition function
    x_new, _ = single(func=acq, method="SLSQP", bounds=bounds, constraints=cons, num_starts=5)

    # evaluate new point
    y_new = func(x_new)

    # add to data
    x_train = torch.vstack((x_train, x_new))
    y_train = torch.hstack((y_train, y_new))

    # print new best
    if y_new > torch.max(y_train[:-1]):
        print(f"New best at evaluation {len(y_train)}: \t Inputs: {x_new.numpy().reshape(dims).round(4)}, \t Outputs: {-y_new.numpy().round(4)}")

# results
best_iter = int(torch.argmax(y_train))
print(f"Evaluation: {best_iter+1} \t Solution: {-float(y_train[best_iter]):.4f}")

New best at evaluation 40:       Inputs: [0.2699 0.     0.2734 0.2569 0.3897 0.5977],    Outputs: [-2.4334]
New best at evaluation 41:       Inputs: [0.2803 0.2137 0.293  0.3056 0.3629 0.5757],    Outputs: [-2.6996]
New best at evaluation 45:       Inputs: [0.2548 0.1987 0.2919 0.2775 0.3478 0.6189],    Outputs: [-2.9417]
New best at evaluation 50:       Inputs: [0.3511 0.1489 0.3001 0.2755 0.3184 0.6503],    Outputs: [-2.9475]
New best at evaluation 51:       Inputs: [0.2964 0.2036 0.4317 0.2867 0.3155 0.642 ],    Outputs: [-3.1666]
New best at evaluation 52:       Inputs: [0.2068 0.1493 0.4381 0.301  0.2881 0.655 ],    Outputs: [-3.2668]
New best at evaluation 54:       Inputs: [0.2156 0.115  0.4715 0.2642 0.3085 0.6714],    Outputs: [-3.2964]
New best at evaluation 57:       Inputs: [0.1959 0.1574 0.4968 0.2744 0.3097 0.6601],    Outputs: [-3.3173]
New best at evaluation 60:       Inputs: [0.2112 0.1458 0.4709 0.2797 0.3103 0.6541],    Outputs: [-3.3201]
New best at evaluation 66:       Inputs: [0.2039 0.142  0.4717 0.2785 0.3116 0.6541],    Outputs: [-3.3209]
Evaluation: 66   Solution: -3.3209