1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124 | import torch
import pandas as pd
import time
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
from warnings import catch_warnings
from warnings import simplefilter
from odyssey.mission import Mission
from odyssey.objective import Objective
from odyssey.navigators import qNParEGO_Navigator
from odyssey.navigators.sampler_navigators import Sobol_Navigator
from botorch.utils.multi_objective.pareto import is_non_dominated
from botorch.utils.multi_objective.box_decompositions.non_dominated import NondominatedPartitioning
def branin(x):
"""
Compute the Branin function
"""
a = 1
b = 5.1 / (4 * torch.pi**2)
c = 5 / torch.pi
r = 6
s = 10
t = 1 / (8 * torch.pi)
x1 = x.squeeze()[0]
x2 = x.squeeze()[1]
return a*(x2 - b*x1**2 + c*x1 - r)**2 + s*(1 - t)*torch.cos(x1) + s
def currin(x):
"""
Compute the Currin function
"""
x1 = x.squeeze()[0]
x2 = x.squeeze()[1]
term1 = 1 - torch.exp(-1 / (2 * x2))
numerator = 2300*x1**3 + 1900*x1**2 + 2092*x1 + 60
denominator = 100*x1**3 + 500*x1**2 + 4*x1 + 20
term2 = numerator / denominator
return term1 + term2
param_space = [(0.0, 1.0), (0.0, 1.0)]
num_iterations = 50
n_init = 15
ref_point = torch.tensor([18.0, 6.0]) # Reference point for hypervolume calculation
############### OPTIMISATION ################
t1 = time.time()
objectives = [
Objective(branin),
Objective(currin),
]
# Set up mission
mission = Mission(
name=f"moo_qnparego",
funcs=objectives,
maneuvers=["descend", "descend"],
envelope=param_space,
)
# Set up navigator
navigator = qNParEGO_Navigator(
mission=mission,
num_init_design=n_init,
input_scaling=False,
data_standardization=False,
init_method=Sobol_Navigator(mission=mission),
acq_function_params={},
)
# BO Loop
opt_samples = len(mission.display_X) - n_init
while opt_samples < num_iterations:
with catch_warnings() as w:
simplefilter('ignore')
trajectory = navigator.trajectory()
observation = navigator.probe(trajectory, init=False)
navigator.relay(trajectory, observation)
navigator.upgrade()
# Compute Pareto front and hypervolume
pareto_mask = is_non_dominated(mission.display_Y)
pareto_Y = mission.display_Y[pareto_mask]
# Update training data
partitioning = NondominatedPartitioning(ref_point=ref_point, Y=pareto_Y)
hypervolume = partitioning.compute_hypervolume().item()
print(f"Iteration {opt_samples+1}: Hypervolume = {hypervolume:.4f}")
opt_samples = len(mission.display_X) - n_init
pareto_Y_final = mission.display_Y[is_non_dominated(mission.display_Y)]
print(f"Finished optimization in {time.time() - t1:.2f} seconds")
################# plot pareto front #################
plt.figure(figsize=(10, 6))
pareto_Y_final_sorted = pareto_Y_final.cpu().numpy()
pareto_Y_ehvi_sorted = pareto_Y_final_sorted[pareto_Y_final_sorted[:, 0].argsort()]
plt.plot(pareto_Y_ehvi_sorted[:, 0], pareto_Y_ehvi_sorted[:, 1], label='qNParEGO Pareto Front', marker='+', markersize=10, color='royalblue', linestyle='-')
plt.scatter(mission.display_Y[:, 0].cpu(), mission.display_Y[:, 1].cpu(), label='qNParEGO Samples', alpha=0.5)
plt.xlabel('Branin')
plt.ylabel('Currin')
plt.title('Pareto Fronts for Branin-Currin using qNParEGO')
plt.legend()
plt.show()
|