# Core Libraries
import numpy as np
import time
import math
import warnings
from scipy.optimize import minimize
from typing import List, Dict, Tuple, Callable, Optional, Union, Any, Iterable
# Scikit-learn for GP, Scaling, and potentially acquisition functions
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, ConstantKernel, WhiteKernel
from sklearn.preprocessing import MinMaxScaler
from sklearn.exceptions import ConvergenceWarning
# SciPy for statistics (needed for Expected Improvement calculation)
from scipy.stats import norm
from functions import bailey_welch_schedule, get_v_star, compute_convolutions, calculate_objective_serv_time_lookup
# Filter warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=ConvergenceWarning) # GP fitting might not always converge perfectly6 Combinatorial Bayesian Optimization Experiments
6.1 Objective
The objective of this experiment is to evaluate and compare the performance of two distinct Combinatorial Bayesian Optimization (CBO) strategies for an outpatient appointment scheduling problem. We investigate:
- CBO utilizing Expected Improvement (EI) as the acquisition function.
- CBO utilizing Lower Confidence Bound (LCB) as the acquisition function with a fixed kappa (\(\kappa\)) value.
We aim to determine which strategy is most effective in identifying an optimal or near-optimal schedule, as measured by the objective function value, leveraging dictionary-based embeddings for the high-dimensional combinatorial space (Deshwal et al. 2023).
6.2 Background
We consider an outpatient appointment scheduling problem as described by Kaandorp and Koole (2007) where the schedule is represented by a vector \(\mathbf{x} = (x_0, x_1, \ldots, x_{T-1})^T\). This vector comprises \(T\) components, where \(x_j\) denotes the non-negative allocation (number of patients) to time slot \(j\), for \(j = 0, \ldots, T-1\). A fundamental constraint is that the total allocation across all time slots must equal a fixed constant \(N\): \[\sum_{j=0}^{T-1} x_j = N\] We require \(x_j \ge 0\) for all \(j = 0, \ldots, T-1\). Consequently, a valid schedule \(\mathbf{x}\) belongs to the feasible set \(\mathcal{F} = \{ \mathbf{z} \in \mathbb{D}^{T} \mid \sum_{j=0}^{T-1} z_j = N, z_j \ge 0 \text{ for all } j\}\), where \(\mathbb{D}\) is the set of non-negative integers (\(\mathbb{Z}_{\ge 0}\)).
Kaandorp and Koole (2007) define a neighborhood structure for local search based on perturbation vectors derived from a set of \(T\) basis change vectors, \(v_i \in \mathbb{D}^{T}\), for \(i = 0, \ldots, T-1\). These basis vectors represent elementary shifts of allocation between time slots:
- \(v_0 = (-1, 0, \ldots, 0, 1)\) (Shift unit from slot 0 to slot \(T-1\))
- \(v_1 = (1, -1, 0, \ldots, 0)\) (Shift unit from slot 1 to slot 0)
- \(v_i = (0, \ldots, 0, \underbrace{1}_{\text{pos } i-1}, \underbrace{-1}_{\text{pos } i}, 0, \ldots, 0)\) for \(i = 2, \ldots, T-1\) (Shift unit from slot \(i\) to slot \(i-1\))
A key property of these basis vectors is that the sum of components for each vector is zero: \(\sum_{j=0}^{T-1} v_{ij} = 0\) for all \(i=0, \ldots, T-1\).
Perturbations are constructed using a binary selection vector \(\mathbf{U} = (u_0, u_1, \ldots, u_{T-1})\), where \(u_i \in \{0, 1\}\). Each \(u_i\) indicates whether the basis change \(v_i\) is included in the perturbation. The resulting perturbation vector \(\mathbf{r}(\mathbf{U}) \in \mathbb{D}^{T}\) is the linear combination: \[\mathbf{r}(\mathbf{U}) := \sum_{i=0}^{T-1} u_i v_i\]
Since each \(v_i\) sums to zero, any perturbation \(\mathbf{r}(\mathbf{U})\) also sums to zero: \(\sum_{j=0}^{T-1} r_j(\mathbf{U}) = 0\). This ensures that applying such a perturbation to a valid schedule \(\mathbf{x}\) preserves the total allocation \(N\).
The neighborhood of a schedule \(\mathbf{x} \in \mathcal{F}\), denoted by \(\mathcal{N}(\mathbf{x})\), comprises all distinct, feasible schedules \(\mathbf{x}'\) reachable by applying a non-zero perturbation \(\mathbf{r}(\mathbf{U})\) (Kaandorp and Koole (2007), use a slightly different but related neighborhood definition based on combinations of these basis vectors).
The objective function to be minimized is a weighted sum of Expected Waiting Time (EWT) and Expected Staff Penalty (ESP), as defined by Kaandorp and Koole (2007): \[C(\mathbf{x}) = w \cdot EWT(\mathbf{x}) + (1-w) \cdot ESP(\mathbf{x})\] Kaandorp and Koole (2007) prove that this objective function is multimodular, which guarantees that a local search algorithm using their defined neighborhood converges to the global optimum.
However, evaluating \(C(\mathbf{x})\) can be computationally expensive, especially for large \(N\) and \(T\). Furthermore, the search space defined by the binary vectors \(\mathbf{U}\) is high-dimensional (\(2^T - 2\) possibilities, excluding \(\mathbf{0}\) and \(\mathbf{1}\)). Bayesian Optimization (BO) is a suitable framework for optimizing such expensive black-box functions. Standard BO methods often struggle with high-dimensional combinatorial spaces. Deshwal et al. (2023) propose a method using dictionary-based embeddings (Hamming Embedding via Dictionaries - HED) to map the high-dimensional binary space of \(\mathbf{U}\) vectors into a lower-dimensional continuous space, where standard Gaussian Process (GP) models can be effectively applied. This experiment applies the HED approach within a BO framework to solve the scheduling problem formulated by Kaandorp and Koole (2007).
6.3 Hypothesis
We hypothesize that:
- Both CBO strategies, leveraging the HED embedding (Deshwal et al. 2023), will be capable of finding schedules superior to the initial schedule derived from the Bailey-Welch method (@).
- CBO strategies employing Lower Confidence Bound (LCB) may exhibit superior performance or faster convergence compared to Expected Improvement (EI), due to the explicit exploration-exploitation trade-off inherent in LCB.
6.4 Methodology
6.4.1 Tools and Materials
- Programming Language: Python 3
- Core Libraries: NumPy, SciPy
- Machine Learning: Scikit-learn (for
GaussianProcessRegressor,MinMaxScaler) - Data Structures: Standard Python lists and dictionaries, NumPy arrays.
- Imported functions:
bailey_welch_schedule,get_v_star,compute_convolutions,calculate_objective_serv_time_lookup(implementing the logic from Bailey (1952), assumed to be in an externalfunctions.pyfile).
6.4.2 Experimental Design
Three distinct Bayesian optimization experiments are conducted, applying the HED embedding approach (Deshwal et al. 2023) to the scheduling problem:
- Experiment 1: Expected Improvement (EI)
- Acquisition Function: Expected Improvement.
- Objective: Minimize \(C(\mathbf{x})\) by iteratively selecting candidate vectors \(\mathbf{U}\) (via their embeddings) that maximize the EI.
- Experiment 2: Lower Confidence Bound (LCB) - Fixed Kappa
- Acquisition Function: Lower Confidence Bound.
- Objective: Minimize \(C(\mathbf{x})\) using a fixed
kappa(\(\kappa\)) value in the LCB acquisition function applied to the GP model over the embedded space.
For all experiments, Hamming Distance Embedding (HED) with a “diverse random” dictionary construction strategy (Deshwal et al. 2023) is employed to map the binary perturbation vectors \(\mathbf{U}\) to a continuous embedding space. A Gaussian Process (GP) model with Automatic Relevance Determination (ARD) kernels models the (negative) objective function in this embedded space.
6.4.3 Variables
- Independent Variables:
- Type of acquisition function (EI, LCB).
- The specific binary perturbation vector \(\mathbf{U}\) selected in each iteration (chosen via optimizing the acquisition function over the embedded space).
- Dependent Variables:
- The objective function value \(C(\mathbf{x}')\) for the resulting schedule \(\mathbf{x}' = \mathbf{x} + \mathbf{r}(\mathbf{U})\) (calculated using the method from Kaandorp and Koole (2007)).
- The best objective function value found throughout the optimization process.
6.4.4 Data Collection
Data, comprising evaluated pairs \((\mathbf{U}, C(\mathbf{x}'))\), is collected iteratively:
- An initial set of
N_INITIALrandomly generated \(\mathbf{U}\) vectors is evaluated. - In each of the subsequent
N_ITERATIONS,BATCH_SIZE_qnew \(\mathbf{U}\) vectors are selected by optimizing the respective acquisition function overNUM_CANDIDATES_Acqfrandomly generated candidate vectors in the original binary space (evaluated via their embeddings). These newly selected vectors are then evaluated, and the results are added to the dataset.
6.4.5 Sample Size and Selection
- N_INITIAL: 20 (number of initial random evaluations)
- N_ITERATIONS: 20 (number of Bayesian optimization iterations)
- BATCH_SIZE_q: 5 (number of candidates selected and evaluated per iteration)
- NUM_CANDIDATES_Acqf: \(T \times 1024 = 20 \times 1024 = 20480\) (number of random candidates generated for optimizing the acquisition function in each iteration)
- m: 128 (dimensionality of the HED embedding space, following Deshwal et al. (2023))
The selection of new points for evaluation is guided by the respective acquisition function (EI or LCB) optimized over the embedded space representation of candidate \(\mathbf{U}\) vectors.
6.4.6 Experimental Procedure
6.4.6.1 1. Setup
Import necessary libraries and configure warning filters.
6.4.6.2 2. Constants
Definition of problem parameters and initial configuration.
# --- Problem Definition ---
# Fixed Data (Use your actual data)
N = 24 # Total number of patients
T = 20 # Dimension of the binary vector U
d = 10 # Length of each interval
max_s = 30 # Maximum service time
q = 0.20 # Probability of a scheduled patient not showing up
w = 0.1 # Weight for the waiting time in objective function
l = 14
v_star = get_v_star(T) # Get the V* matrix(T x T)
# Create service time distribution
def generate_weighted_list(max_s: int, l: float, i: int) -> Optional[np.ndarray]:
"""
Generates a service time probability distribution using optimization.
This function creates a discrete probability distribution over max_s possible
service times (from 1 to max_s). It uses optimization (SLSQP) to find a
distribution whose weighted average service time is as close as possible
to a target value 'l', subject to the constraint that the probabilities
sum to 1 and each probability is between 0 and 1.
After finding the distribution, it sorts the probabilities: the first 'i'
probabilities (corresponding to service times 1 to i) are sorted in
ascending order, and the remaining probabilities (service times i+1 to max_s)
are sorted in descending order.
Note:
- Requires NumPy and SciPy libraries (specifically scipy.optimize.minimize).
Args:
max_s (int): Maximum service time parameter (number of probability bins).
Must be a positive integer.
l (float): The target weighted average service time for the distribution.
Must be between 1 and max_s, inclusive.
i (int): The index determining the sorting split point. Probabilities
for service times 1 to 'i' are sorted ascendingly, and
probabilities for service times 'i+1' to 'max_s' are sorted
descendingly. Must be between 1 and max_s-1 for meaningful sorting.
Returns:
numpy.ndarray: An array of size max_s+1. The first element (index 0) is 0.
Elements from index 1 to max_s represent the calculated
and sorted probability distribution, summing to 1.
Returns None if optimization fails or inputs are invalid.
"""
# --- Input Validation ---
if not isinstance(max_s, int) or max_s <= 0:
print(f"Error: max_s must be a positive integer, but got {max_s}")
return None
if not isinstance(l, (int, float)) or not (1 <= l <= max_s):
print(f"Error: Target average 'l' ({l}) must be between 1 and max_s ({max_s}).")
return None
if not isinstance(i, int) or not (0 < i < max_s):
print(f"Error: Sorting index 'i' ({i}) must be between 1 and max_s-1 ({max_s-1}).")
# If clamping is desired instead of error:
# print(f"Warning: Index 'i' ({i}) is outside the valid range (1 to {max_s-1}). Clamping i.")
# i = max(1, min(i, max_s - 1))
return None # Strict check based on docstring requirement
# --- Inner helper function for optimization ---
def objective(x: np.ndarray) -> float:
"""Objective function: Squared difference between weighted average and target l."""
# x represents probabilities P(1) to P(max_s)
service_times = np.arange(1, max_s + 1)
weighted_avg = np.dot(service_times, x) # Equivalent to sum(k * P(k) for k=1 to max_s)
return (weighted_avg - l) ** 2
# --- Constraints for optimization ---
# Constraint 1: The sum of the probabilities must be 1
constraints = ({
'type': 'eq',
'fun': lambda x: np.sum(x) - 1.0 # Ensure float comparison
})
# Bounds: Each probability value x[k] must be between 0 and 1
# Creates a list of max_s tuples, e.g., [(0, 1), (0, 1), ..., (0, 1)]
bounds = [(0, 1)] * max_s
# Initial guess: Use Dirichlet distribution to get a random distribution that sums to 1.
# Provides a starting point for the optimizer. np.ones(max_s) gives equal weights initially.
initial_guess = np.random.dirichlet(np.ones(max_s))
# --- Perform Optimization ---
try:
result = minimize(
objective,
initial_guess,
method='SLSQP',
bounds=bounds,
constraints=constraints,
# options={'disp': False} # Set True for detailed optimizer output
)
# Check if optimization was successful
if not result.success:
print(f"Warning: Optimization failed! Message: {result.message}")
# Optionally print result object for more details: print(result)
return None # Indicate failure
# The optimized probabilities (P(1) to P(max_s))
optimized_probs = result.x
# --- Post-process: Correct potential floating point inaccuracies ---
# Ensure probabilities are non-negative and sum *exactly* to 1
optimized_probs[optimized_probs < 0] = 0 # Clamp small negatives to 0
current_sum = np.sum(optimized_probs)
if not np.isclose(current_sum, 1.0):
if current_sum > 0: # Avoid division by zero
optimized_probs /= current_sum # Normalize to sum to 1
else:
print("Warning: Optimization resulted in zero sum probabilities after clamping negatives.")
# Handle this case - maybe return uniform distribution or None
return None # Or return uniform: np.ones(max_s) / max_s
except Exception as e:
print(f"An error occurred during optimization: {e}")
return None
# --- Reorder the probabilities based on the index 'i' ---
# Split the probabilities P(1)...P(i) and P(i+1)...P(max_s)
# Note: Python slicing is exclusive of the end index, array indexing is 0-based.
# result.x[0] corresponds to P(1), result.x[i-1] to P(i).
# result.x[i] corresponds to P(i+1), result.x[max_s-1] to P(max_s).
first_part_probs = optimized_probs[:i] # Probabilities P(1) to P(i)
second_part_probs = optimized_probs[i:] # Probabilities P(i+1) to P(max_s)
# Sort the first part ascending, the second part descending
sorted_first_part = np.sort(first_part_probs)
sorted_second_part = np.sort(second_part_probs)[::-1] # [::-1] reverses
# --- Create final output array ---
# Array of size max_s + 1, initialized to zeros. Index 0 unused.
values = np.zeros(max_s + 1)
# Assign the sorted probabilities back into the correct slots (index 1 onwards)
values[1 : i + 1] = sorted_first_part # Assign P(1)...P(i)
values[i + 1 : max_s + 1] = sorted_second_part # Assign P(i+1)...P(max_s)
# Final check on sum after potential normalization/sorting
if not np.isclose(np.sum(values[1:]), 1.0):
print(f"Warning: Final distribution sum is {np.sum(values[1:])}, not 1.0. Check logic.")
# Return the final array with the sorted probability distribution
return values
i = 10 # First 5 highest values in ascending order, rest in descending order
s = generate_weighted_list(max_s, l, i)
print(f"Average generated service time: {np.dot(np.arange(len(s)), s)}")
convolutions = compute_convolutions(s, N, q)
X = np.array(bailey_welch_schedule(T, d, N, s))
print(f"Initial schedule: {X}")
# Objective Function Calculation
LARGE_PENALTY = 1e10 # Penalty for infeasible solutions
ewt, esp = calculate_objective_serv_time_lookup(X, d, convolutions)
initial_objective_value = w * ewt + (1 - w) * esp
print(f"Initial objective value: {initial_objective_value}")Average generated service time: 12.942391896136673
Initial schedule: [2 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1 1 0 1 8]
Initial objective value: 120.674268580054476.4.6.3 3. Common Functions (Objective Evaluation and HED)
Objective evaluation implements \(C(\mathbf{x})\) from Kaandorp and Koole (2007). HED implementation follows Deshwal et al. (2023).
def evaluate_objective(U_np, X_vec, v_star, convolutions, d, w):
"""
Target function: Evaluates objective for a single binary numpy array U.
Returns a float.
"""
# Input validation (same as before)
if not isinstance(U_np, np.ndarray):
raise TypeError("Input U must be a numpy array")
if U_np.ndim != 1:
raise ValueError("Input U must be 1-dimensional")
if U_np.shape[0] != v_star.shape[0]:
raise ValueError(f"Dimension mismatch: U length {U_np.shape[0]} != V* rows {v_star.shape[0]}.")
if X_vec.shape[0] != v_star.shape[1]:
raise ValueError("Dimension mismatch: X length must match V* columns.")
if not np.all((U_np == 0) | (U_np == 1)):
raise ValueError("Input U must be binary (0s and 1s).")
# Calculate Y based on selected rows of V_star
V_sum = np.sum(v_star[U_np == 1, :], axis=0)
Y = X_vec + V_sum
# Check feasibility and calculate objective
if np.all(Y >= 0):
ewt, esp = calculate_objective_serv_time_lookup(Y, d, convolutions)
objective_value = w * ewt + (1 - w) * esp
return objective_value
else:
# Infeasible solution
return LARGE_PENALTY
# --- HED Implementation ---
def hamming_distance(u1, u2):
"""Calculates Hamming distance between two binary numpy arrays."""
return np.sum(u1 != u2)
def generate_diverse_random_dictionary(T, m):
"""Generates the random dictionary A for HED."""
dictionary_A = np.zeros((m, T), dtype=int)
for i in range(m):
# Sample theta for density of 1s in this dictionary vector
theta = np.random.uniform(0, 1)
row = (np.random.rand(T) < theta).astype(int)
dictionary_A[i, :] = row
return dictionary_A
def _generate_binary_hadamard_matrix_recursive(dim):
"""
Generates a binary (0/1) Hadamard-like matrix of size dim x dim.
'dim' must be a power of 2.
This uses the Sylvester's construction H_2n = [[H_n, H_n], [H_n, 1-H_n]]
starting with H_1 = [[1]].
"""
if not (dim > 0 and (dim & (dim - 1) == 0)): # Checks if dim is a power of 2
raise ValueError("Dimension must be a power of 2.")
if dim == 1:
return np.array([[1]], dtype=int)
else:
h_prev = _generate_binary_hadamard_matrix_recursive(dim // 2)
h_top = np.hstack((h_prev, h_prev))
h_bottom = np.hstack((h_prev, 1 - h_prev)) # 1-H_n for binary
return np.vstack((h_top, h_bottom))
def generate_wavelet_dictionary(T, m):
"""
Generates a dictionary A of size m x T using the subsampled binary wavelet approach.
Args:
T (int): The dimensionality of the input space (number of columns in dictionary).
m (int): The desired number of dictionary elements (number of rows).
Returns:
np.ndarray: An m x T integer numpy array representing the dictionary.
"""
if T <= 0:
raise ValueError("T (dimensionality) must be positive.")
if m <= 0:
raise ValueError("m (dictionary size) must be positive.")
# 1. Determine the smallest power of 2 >= T for the full wavelet matrix
if T == 1:
n_wavelet = 1
elif (T > 0 and (T & (T - 1) == 0)): # T is already a power of 2
n_wavelet = T
else:
n_wavelet = 2**math.ceil(math.log2(T))
# 2. Generate the full n_wavelet x n_wavelet binary Hadamard matrix
# print(f"Generating full wavelet matrix of size: {n_wavelet}x{n_wavelet}")
full_wavelet_matrix = _generate_binary_hadamard_matrix_recursive(n_wavelet)
# 3. Subsample T columns if n_wavelet > T
if n_wavelet > T:
# print(f"Subsampling {T} columns from {n_wavelet} columns.")
col_indices = np.random.choice(n_wavelet, size=T, replace=False)
col_indices.sort() # Optional: for deterministic testing if seed is set
wavelet_matrix_T_cols = full_wavelet_matrix[:, col_indices]
else:
wavelet_matrix_T_cols = full_wavelet_matrix # n_wavelet == T
# 4. Subsample m rows
num_available_rows = wavelet_matrix_T_cols.shape[0]
if m > num_available_rows:
print(f"Warning: Requested dictionary size m ({m}) is greater than "
f"available unique wavelet rows ({num_available_rows}). "
f"Using all available rows and repeating if necessary, or consider reducing m.")
# For simplicity, if m > num_available_rows, we'll sample with replacement
# or you could choose to error, or return fewer rows.
# The paper implies m should be less than or equal to the rows of B_d.
# If sampling with replacement is needed:
row_indices = np.random.choice(num_available_rows, size=m, replace=True)
# If strictly no replacement and m > num_available_rows, one might error or cap m
# row_indices = np.random.choice(num_available_rows, size=min(m, num_available_rows), replace=False)
# if m > num_available_rows:
# # Handle the case where more rows are needed than available unique ones
# # This might involve repeating rows or another strategy
# pass
else:
# print(f"Subsampling {m} rows from {num_available_rows} available rows.")
row_indices = np.random.choice(num_available_rows, size=m, replace=False)
row_indices.sort() # Optional: for deterministic testing if seed is set
dictionary_A = wavelet_matrix_T_cols[row_indices, :]
return dictionary_A
def embed_vector(U_np, dictionary_A):
"""Embeds a single binary vector U using HED."""
m = dictionary_A.shape[0]
embedding_phi = np.zeros(m, dtype=float) # Use float for GP
for i in range(m):
embedding_phi[i] = hamming_distance(U_np, dictionary_A[i, :])
return embedding_phi
def embed_batch(U_batch_np, dictionary_A):
"""Embeds a batch of binary vectors U."""
# Input U_batch_np is expected to be a NumPy array
m = dictionary_A.shape[0]
if U_batch_np.ndim == 1: # Handle single vector case
U_batch_np = U_batch_np.reshape(1, -1)
batch_size = U_batch_np.shape[0]
embeddings_np = np.zeros((batch_size, m), dtype=float) # Use float for GP
for j in range(batch_size):
embeddings_np[j, :] = embed_vector(U_batch_np[j, :], dictionary_A)
# Return NumPy array directly
return embeddings_np6.4.6.4 4. Experiment 1: CBO with Expected Improvement (EI)
Applies the methodology from Deshwal et al. (2023) using EI.
# --- BO Helper Functions ---
def get_fitted_model(train_X_embedded_scaled, train_Y, m):
"""
Fits a GaussianProcessRegressor model to the SCALED embedded data.
Assumes train_Y contains negative objective values for maximization.
"""
if train_Y.ndim > 1 and train_Y.shape[1] == 1:
train_Y = train_Y.ravel() # sklearn GP expects 1D target array
# Define the kernel for the Gaussian Process
# Matern kernel is a common choice, nu=2.5 is smooth (twice differentiable)
# ConstantKernel handles the overall variance scaling
# WhiteKernel handles the observation noise
kernel = ConstantKernel(1.0, constant_value_bounds=(1e-3, 1e3)) * \
Matern(length_scale=np.ones(m), # Enable ARD, initialize length scales to 1
length_scale_bounds=(1e-2, 1e2),
nu=2.5) + \
WhiteKernel(noise_level=1e-10, # Small value for numerical stability
noise_level_bounds="fixed") # Bounds for noise optimization
# Instantiate the Gaussian Process Regressor
# alpha: Value added to the diagonal of the kernel matrix during fitting
# for numerical stability (can also be seen as additional noise)
# n_restarts_optimizer: Restarts optimizer to find better hyperparameters
gp_model = GaussianProcessRegressor(
kernel=kernel,
alpha=1e-10, # Small value for numerical stability
n_restarts_optimizer=10, # More restarts -> better hyperparams but slower
random_state=42 # For reproducibility of optimizer restarts
)
# Fit the GP model
gp_model.fit(train_X_embedded_scaled, train_Y)
return gp_model
def expected_improvement(mu, sigma, f_best, xi=0.01):
"""
Computes the Expected Improvement acquisition function.
Assumes maximization (f_best is the current maximum observed value).
mu, sigma: Predicted mean and standard deviation (NumPy arrays).
f_best: Current best observed function value (scalar).
xi: Exploration-exploitation trade-off parameter.
"""
# Ensure sigma is positive and non-zero to avoid division errors
sigma = np.maximum(sigma, 1e-9)
Z = (mu - f_best - xi) / sigma
ei = (mu - f_best - xi) * norm.cdf(Z) + sigma * norm.pdf(Z)
# Set EI to 0 where variance is negligible
ei[sigma <= 1e-9] = 0.0
return ei
# MODIFIED: Accepts the scaler and uses scikit-learn GP + EI
def optimize_acqf_discrete_via_embedding(gp_model, scaler, dictionary_A, T, q, num_candidates, current_best_neg_f_val):
"""
Optimizes acquisition function (Expected Improvement) by sampling random
binary candidates, embedding, SCALING, predicting with GP, and calculating EI.
Selects the top q candidates based on EI.
Returns candidates as a numpy array (q x T).
"""
m = dictionary_A.shape[0]
# 1. Generate Random Binary Candidates
candidate_u_vectors_np = np.random.randint(0, 2, size=(num_candidates, T))
# Optional: Ensure unique candidates if needed (adds overhead)
# candidate_u_vectors_np = np.unique(candidate_u_vectors_np, axis=0)
# num_candidates = candidate_u_vectors_np.shape[0] # Update count
# 2. Embed the Candidates
embedded_candidates_np = embed_batch(candidate_u_vectors_np, dictionary_A)
# 3. Scale the Embedded Candidates
# Handle potential warning if scaler expects float64 (already float here)
# Use the *fitted* scaler from the training data
embedded_candidates_scaled_np = scaler.transform(embedded_candidates_np)
# 4. Predict Mean and Std Dev using the GP Model
mu, std = gp_model.predict(embedded_candidates_scaled_np, return_std=True)
# 5. Calculate Acquisition Function (Expected Improvement)
# current_best_neg_f_val is the maximum of the (negative) objectives seen so far
acq_values = expected_improvement(mu, std, current_best_neg_f_val, xi=0.01)
# 6. Select Top Candidates
# Use np.argsort to find indices that would sort the array (ascending)
# Select the last q indices for the highest EI values
# If q=1, np.argmax(acq_values) is simpler but argsort works generally
top_indices = np.argsort(acq_values)[-q:]
# Ensure indices are returned in descending order of acquisition value (optional but nice)
top_indices = top_indices[::-1]
return candidate_u_vectors_np[top_indices, :]
# --- BO Loop ---
# Parameters
N_INITIAL = 49
N_ITERATIONS = 25
BATCH_SIZE_q = 5
NUM_CANDIDATES_Acqf = T*3*1024 # Might need more for higher T
m = math.ceil(T/4) # Dimension of the embedding space
# Store evaluated points (using NumPy arrays)
evaluated_U_np_list = [] # List to store evaluated U vectors (binary)
evaluated_f_vals = [] # List to store raw objective values (lower is better)
train_Y_list = [] # List to store NEGATED objective values for GP (higher is better)
# 1. Initialization
print(f"Generating {N_INITIAL} initial points...")
initial_candidates = [np.zeros(T, dtype=int)]
while len(initial_candidates) < N_INITIAL + 1: # +1 for the zero vector = initial X
U_init = np.random.randint(0, 2, size=T)
# Ensure unique initial points
is_duplicate = any(np.array_equal(U_init, u) for u in initial_candidates)
if not is_duplicate:
initial_candidates.append(U_init)
for U_init in initial_candidates:
f_val = evaluate_objective(U_init, X, v_star, convolutions, d, w)
neg_f_val = -f_val
evaluated_U_np_list.append(U_init)
evaluated_f_vals.append(f_val)
train_Y_list.append(neg_f_val)
# Convert lists to NumPy arrays for GP fitting
train_Y = np.array(train_Y_list).reshape(-1, 1) # Keep as column vector initially
best_obj_so_far = min(evaluated_f_vals) if evaluated_f_vals else float('inf')
initial_best_obj_so_far_ei = best_obj_so_far
print(f"Initial best objective value: {best_obj_so_far}")
if not np.isfinite(best_obj_so_far):
print("Warning: Initial best objective is infinite, possibly all initial points were infeasible.")
# 2. BO Iterations
for iteration in range(N_ITERATIONS):
start_time = time.time()
print(f"\n--- Iteration {iteration + 1}/{N_ITERATIONS} ---")
# a. Generate dictionary A for HED
current_dictionary_A = generate_diverse_random_dictionary(T, m)
# b. Embed ALL evaluated U vectors so far
if not evaluated_U_np_list:
print("Warning: No points evaluated yet. Skipping iteration.")
continue
evaluated_U_np_array = np.array(evaluated_U_np_list)
embedded_train_X = embed_batch(evaluated_U_np_array, current_dictionary_A)
# c. Scale the embedded training data
scaler = MinMaxScaler()
# Fit scaler only if there's data
if embedded_train_X.shape[0] > 0:
# Fit and transform
embedded_train_X_scaled = scaler.fit_transform(embedded_train_X)
else:
# Handle case with no data (shouldn't happen after init)
embedded_train_X_scaled = embedded_train_X # Will be empty
# Ensure train_Y is a NumPy array for fitting
train_Y_for_fit = np.array(train_Y_list) # Use the list directly
# d. Fit GP Model using SCALED data
print("Fitting GP model...")
if embedded_train_X_scaled.shape[0] > 0 and train_Y_for_fit.shape[0] == embedded_train_X_scaled.shape[0]:
gp_model = get_fitted_model(embedded_train_X_scaled, train_Y_for_fit, m)
print("GP model fitted.")
else:
print("Warning: Not enough data or data mismatch to fit GP model. Skipping iteration.")
continue # Skip if no data or mismatch
# e. Determine current best value for Acquisition Function
# We are maximizing the negative objective in the GP
current_best_neg_f_val = np.max(train_Y_for_fit) if train_Y_for_fit.size > 0 else -float('inf')
# Prevent potential issues if all points were infeasible (very large negative best_f)
if current_best_neg_f_val <= -LARGE_PENALTY / 2 and np.isfinite(current_best_neg_f_val):
print(f"Warning: Current best value ({current_best_neg_f_val:.2f}) is very low (likely from penalties). Acqf might behave unexpectedly.")
# f. Optimize Acquisition Function (Expected Improvement)
print("Optimizing acquisition function...")
next_U_candidates_np = optimize_acqf_discrete_via_embedding(
gp_model=gp_model,
scaler=scaler, # Pass the fitted scaler
dictionary_A=current_dictionary_A,
T=T,
q=BATCH_SIZE_q,
num_candidates=NUM_CANDIDATES_Acqf,
current_best_neg_f_val=current_best_neg_f_val
)
print(f"Selected {next_U_candidates_np.shape[0]} candidate(s).")
# g. Evaluate Objective for the selected candidate(s)
newly_evaluated_U = []
newly_evaluated_f = []
newly_evaluated_neg_f = []
for i in range(next_U_candidates_np.shape[0]):
next_U = next_U_candidates_np[i, :]
# Check if this candidate was already evaluated
# Use a tolerance for floating point comparisons if U were continuous
# For binary, exact comparison is fine
already_evaluated = any(np.array_equal(next_U, u) for u in evaluated_U_np_list)
if already_evaluated:
print(f" Candidate {i} was already evaluated. Skipping re-evaluation.")
# TODO: Optionally, could try to generate a *different* candidate here
# e.g., by running optimize_acqf again excluding this one,
# or sampling randomly near it. For now, just skip.
continue # Skip to next candidate
# Evaluate the objective
next_f = evaluate_objective(next_U, X, v_star, convolutions, d, w)
next_neg_f = -next_f
V_sum = np.sum(v_star[next_U == 1, :], axis=0)
Y = X + V_sum
print(f" Candidate {i}: Obj = {next_f:.4f}, schedule = {Y}")
# Add to temporary lists for this iteration
newly_evaluated_U.append(next_U)
newly_evaluated_f.append(next_f)
newly_evaluated_neg_f.append(next_neg_f)
# Update overall best objective found
if next_f < best_obj_so_far:
best_obj_so_far = next_f
# h. Augment Dataset for next iteration
evaluated_U_np_list.extend(newly_evaluated_U)
evaluated_f_vals.extend(newly_evaluated_f)
train_Y_list.extend(newly_evaluated_neg_f) # Add negative values for next GP fit
# Convert train_Y_list back to array for potential use (though we rebuild it next iter)
train_Y = np.array(train_Y_list).reshape(-1, 1)
iter_time = time.time() - start_time
print(f"Best objective value found so far: {best_obj_so_far:.4f}")
print(f"Total points evaluated: {len(evaluated_f_vals)}")
print(f"Iteration {iteration + 1} completed in {iter_time:.2f} seconds.")
# --- Results ---
print("\n--- Optimization Finished ---")
if not evaluated_f_vals:
print("No points were successfully evaluated.")
else:
# Find the best point among all evaluated points
final_best_idx_ei = np.argmin(evaluated_f_vals) # Index of minimum raw objective
final_best_U_ei = evaluated_U_np_list[final_best_idx_ei]
final_best_f_ei = evaluated_f_vals[final_best_idx_ei]
nr_evaluated_f_vals_ei = len(evaluated_f_vals) # Saved for reporting results
print(f"Total evaluations: {len(evaluated_f_vals)}")
print(f"Best Objective Value Found: {final_best_f_ei}")
# Ensure U is printed correctly if it's long
print(f"Best U vector Found: {final_best_U_ei}")
# print(f"Best U vector Found (Indices of 1s): {np.where(final_best_U_ei == 1)[0]}")
# Verification - Recalculate Y for the best U found
V_sum_best = np.sum(v_star[final_best_U_ei == 1, :], axis=0)
Y_best_ei = X + V_sum_best
is_feasible = np.all(Y_best_ei >= 0)
if is_feasible:
ewt, esp = calculate_objective_serv_time_lookup(Y_best_ei, d, convolutions)
recalculated_obj = w * ewt + (1 - w) * esp
else:
LARGE_PENALTY
print(f"\n--- Verification ---")
print(f"Is the best U feasible? {is_feasible}")
if is_feasible:
print(f"Resulting Y vector for best U: {Y_best_ei}")
print(f"Objective value (recalculated): {recalculated_obj:.4f}")
if not np.isclose(final_best_f_ei, recalculated_obj):
print(f"Warning: Stored best objective ({final_best_f_ei}) does not match recalculation ({recalculated_obj})!")
elif final_best_f < LARGE_PENALTY:
print(f"Warning: Best objective ({final_best_f_ei}) is not the penalty value, but feasibility check failed.")
print(f"Resulting Y vector (infeasible): {Y_best_ei}")
else:
print("Best solution found corresponds to an infeasible penalty value.")Generating 49 initial points...
Initial best objective value: 118.87100166955729
--- Iteration 1/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 114.9636, schedule = [2 1 1 1 1 1 1 0 1 1 1 1 0 1 1 0 2 0 1 7]
Candidate 1: Obj = 130.2052, schedule = [1 1 1 1 1 1 1 0 1 1 0 2 1 0 1 0 2 0 0 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 2 0 1 1 1 1 -1 2 0 2 1 0 1 1 0 1 0 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 1 1 1 0 1 2 1 0 1 0 2 1 0 1 0 2 -1 1 9]
Candidate 4: Obj = 131.1795, schedule = [1 1 2 1 0 1 1 1 0 1 0 2 0 1 1 0 2 0 0 9]
Best objective value found so far: 114.9636
Total points evaluated: 55
Iteration 1 completed in 2.91 seconds.
--- Iteration 2/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 124.1641, schedule = [3 1 1 1 0 1 1 0 1 0 1 1 1 0 1 1 2 0 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 0 1 2 0 1 1 1 -1 1 1 1 2 0 0 1 2 0 0 8]
Candidate 2: Obj = 123.1027, schedule = [3 0 1 2 0 1 1 1 0 0 1 1 1 1 0 1 2 0 0 8]
Candidate 3: Obj = 124.6251, schedule = [3 0 2 1 0 1 0 1 0 1 1 2 1 0 0 1 2 0 0 8]
Candidate 4: Obj = 124.1919, schedule = [3 1 0 2 0 1 1 0 1 0 2 0 1 0 1 1 2 0 0 8]
Best objective value found so far: 114.9636
Total points evaluated: 60
Iteration 2 completed in 2.75 seconds.
--- Iteration 3/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 0 1 1 1 1 -1 2 1 0 2 -1 1 2 1 -1 2 8]
Candidate 1: Obj = 126.5321, schedule = [2 1 1 0 1 1 0 2 0 1 1 0 1 0 1 2 1 0 1 8]
Candidate 2: Obj = 122.9236, schedule = [2 1 2 0 1 1 1 0 1 0 2 1 1 0 1 0 2 0 0 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 1 1 0 0 2 1 0 1 0 1 2 1 -1 2 0 2 0 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 1 0 0 2 0 2 1 -1 1 1 2 0 1 7]
Best objective value found so far: 114.9636
Total points evaluated: 65
Iteration 3 completed in 2.41 seconds.
--- Iteration 4/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 2 1 -1 2 0 1 0 2 0 2 1 0 0 1 2 -1 2 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 2 1 0 1 1 0 0 2 1 1 1 0 0 1 2 -1 1 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 0 1 1 1 1 1 1 1 0 0 1 2 -1 2 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 0 1 1 0 2 0 1 1 1 1 1 1 0 0 1 2 -1 2 8]
Candidate 4: Obj = 132.3160, schedule = [2 0 2 1 0 1 0 1 1 1 0 2 1 0 0 1 2 0 0 9]
Best objective value found so far: 114.9636
Total points evaluated: 70
Iteration 4 completed in 2.65 seconds.
--- Iteration 5/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 1 1 0 1 1 1 1 -1 2 1 1 0 1 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 2 0 1 1 1 0 1 2 -1 1 1 1 1 1 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 1 1 1 1 0 2 0 2 -1 2 0 2 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 1 1 1 0 0 2 1 -1 1 1 2 0 1 0 2 0 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 1 1 0 1 2 -1 2 1 0 1 0 1 2 1 -1 2 7]
Best objective value found so far: 114.9636
Total points evaluated: 75
Iteration 5 completed in 2.52 seconds.
--- Iteration 6/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 2 0 0 2 0 2 -1 2 1 1 1 0 0 1 1 1 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 2 1 -1 1 2 0 1 1 0 1 1 1 0 1 1 1 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 1 2 0 0 2 0 2 -1 2 1 1 1 0 1 1 0 0 1 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 2 1 1 0 1 0 1 1 1 0 1 2 0 0 2 1 -1 2 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 1 1 0 1 1 1 1 1 0 2 -1 2 0 1 1 0 8]
Best objective value found so far: 114.9636
Total points evaluated: 80
Iteration 6 completed in 2.68 seconds.
--- Iteration 7/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 116.6354, schedule = [3 0 2 0 1 1 0 1 1 1 1 1 1 0 1 0 2 0 1 7]
Candidate 1: Obj = 116.6189, schedule = [3 1 1 1 0 1 0 1 1 1 1 0 2 0 1 1 1 0 1 7]
Candidate 2: Obj = 124.8705, schedule = [2 0 1 2 0 1 1 0 1 1 1 1 0 1 0 1 2 0 1 8]
Candidate 3: Obj = 123.0720, schedule = [3 1 1 1 0 1 0 1 1 1 1 1 1 0 1 1 1 0 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 1 1 1 1 1 0 1 1 1 1 -1 2 0 2 -1 2 8]
Best objective value found so far: 114.9636
Total points evaluated: 85
Iteration 7 completed in 2.61 seconds.
--- Iteration 8/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 125.8226, schedule = [2 1 1 0 1 0 1 2 0 0 2 0 1 1 1 0 2 0 1 8]
Candidate 1: Obj = 122.2661, schedule = [2 1 1 1 1 1 0 1 0 1 2 0 1 0 2 1 1 0 0 8]
Candidate 2: Obj = 126.1646, schedule = [3 1 1 0 0 2 0 1 0 1 1 2 1 0 0 2 0 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 1 0 2 0 1 0 2 1 -1 2 0 1 1 0 9]
Candidate 4: Obj = 125.4905, schedule = [2 1 2 1 0 0 1 1 0 2 0 2 1 0 0 2 0 0 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 90
Iteration 8 completed in 2.83 seconds.
--- Iteration 9/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 0 2 1 -1 2 0 1 1 0 2 0 2 -1 2 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 1 2 1 0 0 2 1 -1 2 0 1 1 0 2 0 2 -1 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 0 1 2 0 0 2 0 0 2 0 2 1 -1 1 2 1 -1 1 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 0 1 2 0 0 2 1 -1 1 1 2 0 0 1 2 1 -1 2 8]
Candidate 4: Obj = 132.6410, schedule = [2 0 2 1 0 0 2 0 1 1 0 2 0 0 1 1 1 0 2 8]
Best objective value found so far: 114.9636
Total points evaluated: 95
Iteration 9 completed in 2.84 seconds.
--- Iteration 10/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 122.2296, schedule = [2 2 0 1 1 0 2 0 1 1 0 2 0 1 0 2 0 0 2 7]
Candidate 1: Obj = 119.2311, schedule = [3 1 0 1 1 0 2 0 0 1 1 2 1 0 0 1 1 1 1 7]
Candidate 2: Obj = 121.8186, schedule = [3 1 0 1 0 1 1 1 1 0 1 2 1 0 1 1 0 1 0 8]
Candidate 3: Obj = 130.3095, schedule = [1 2 0 1 1 1 0 1 1 0 1 2 0 1 1 0 1 0 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 2 0 0 1 1 1 0 1 2 1 -1 1 2 0 1 1 7]
Best objective value found so far: 114.9636
Total points evaluated: 100
Iteration 10 completed in 2.78 seconds.
--- Iteration 11/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 1 2 0 1 1 0 1 0 2 0 2 -1 1 1 2 -1 1 9]
Candidate 1: Obj = 132.3651, schedule = [2 0 1 1 1 1 1 0 0 2 0 2 0 1 0 1 2 0 0 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 0 1 1 1 1 0 1 1 1 0 2 1 -1 2 0 2 -1 1 9]
Candidate 3: Obj = 132.7032, schedule = [1 1 1 2 0 0 2 0 1 1 0 2 0 0 2 0 2 0 0 9]
Candidate 4: Obj = 133.3754, schedule = [2 0 1 2 0 1 1 0 0 2 0 2 0 0 2 0 1 1 0 9]
Best objective value found so far: 114.9636
Total points evaluated: 105
Iteration 11 completed in 2.70 seconds.
--- Iteration 12/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 1 -1 2 1 1 0 0 1 1 1 1 1 0 2 -1 1 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 1 0 0 2 0 2 0 0 1 2 1 0 1 0 2 -1 1 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 0 2 0 1 0 1 2 -1 1 2 1 0 1 0 2 1 0 1 8]
Candidate 3: Obj = 127.0694, schedule = [2 0 2 0 0 1 2 1 0 0 1 2 0 1 1 0 2 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 2 1 1 0 1 1 1 0 0 1 1 1 1 0 1 2 -1 2 8]
Best objective value found so far: 114.9636
Total points evaluated: 110
Iteration 12 completed in 2.93 seconds.
--- Iteration 13/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 134.6761, schedule = [1 2 1 0 0 1 2 0 0 1 1 2 0 0 1 1 1 0 1 9]
Candidate 1: Obj = 134.2275, schedule = [1 1 1 2 0 0 2 0 0 2 0 1 1 0 2 0 1 0 1 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 2 0 0 2 1 1 -1 1 1 2 1 -1 1 1 1 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 2 0 1 0 2 0 2 -1 1 1 1 1 0 2 0 1 1 0 9]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 1 2 1 -1 2 0 1 1 0 1 2 0 1 0 9]
Best objective value found so far: 114.9636
Total points evaluated: 115
Iteration 13 completed in 3.68 seconds.
--- Iteration 14/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 2 1 1 0 0 1 2 -1 1 2 1 0 0 2 0 1 0 2 7]
Candidate 1: Obj = 122.5544, schedule = [3 1 1 0 1 0 1 1 0 2 0 1 1 0 1 1 1 1 0 8]
Candidate 2: Obj = 119.2175, schedule = [2 1 2 0 1 1 1 0 0 2 1 0 1 0 1 1 2 0 1 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 2 -1 2 0 1 1 0 2 0 2 -1 2 0 1 1 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 2 1 0 1 0 2 0 1 1 1 0 2 -1 1 1 1 1 1 7]
Best objective value found so far: 114.9636
Total points evaluated: 120
Iteration 14 completed in 3.42 seconds.
--- Iteration 15/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 2 0 1 1 1 1 -1 2 1 1 0 0 2 0 2 0 0 9]
Candidate 1: Obj = 133.5979, schedule = [2 0 2 0 0 2 1 0 0 1 1 2 1 0 1 0 1 0 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 2 1 1 0 1 0 1 1 1 1 0 1 0 2 1 1 -1 2 7]
Candidate 3: Obj = 125.6929, schedule = [3 1 1 0 1 1 1 0 1 0 2 1 1 0 1 0 1 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 1 1 1 0 1 1 1 1 0 2 -1 2 1 0 1 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 125
Iteration 15 completed in 4.66 seconds.
--- Iteration 16/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 131.6450, schedule = [2 0 1 1 0 1 2 0 1 1 0 1 2 0 0 1 2 0 0 9]
Candidate 1: Obj = 131.8712, schedule = [2 1 0 1 0 2 0 1 1 1 0 1 2 0 0 1 2 0 0 9]
Candidate 2: Obj = 132.0733, schedule = [2 1 0 1 0 2 0 2 0 1 0 1 2 0 1 0 2 0 0 9]
Candidate 3: Obj = 130.7426, schedule = [2 0 1 1 0 1 2 0 1 1 0 1 2 0 1 1 1 0 0 9]
Candidate 4: Obj = 131.2668, schedule = [2 1 0 1 0 1 1 2 0 1 0 1 2 0 1 1 1 0 0 9]
Best objective value found so far: 114.9636
Total points evaluated: 130
Iteration 16 completed in 3.87 seconds.
--- Iteration 17/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 2 1 -1 2 1 1 -1 1 1 2 0 0 1 1 2 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 1 2 0 1 0 2 -1 1 1 2 0 1 0 2 1 -1 1 9]
Candidate 2: Obj = 126.8993, schedule = [1 1 1 1 1 1 0 1 0 2 0 2 1 0 0 2 0 1 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 1 1 1 1 1 -1 1 2 1 0 0 1 2 1 0 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 2 0 0 2 0 2 -1 2 0 2 1 -1 1 1 2 0 1 7]
Best objective value found so far: 114.9636
Total points evaluated: 135
Iteration 17 completed in 3.84 seconds.
--- Iteration 18/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 2 0 0 1 2 0 0 1 2 0 1 1 1 0 2 -1 2 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 0 2 0 0 1 2 0 0 2 1 1 0 1 0 1 2 -1 2 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 2 0 0 1 2 1 -1 1 2 1 0 1 1 0 2 -1 2 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 2 0 1 1 0 2 1 1 0 1 0 2 -1 2 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 2 0 0 1 2 1 -1 1 2 0 1 1 1 0 2 -1 2 7]
Best objective value found so far: 114.9636
Total points evaluated: 140
Iteration 18 completed in 4.55 seconds.
--- Iteration 19/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 1 1 2 -1 2 1 0 1 0 1 1 2 -1 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 0 1 2 0 1 1 1 0 0 1 1 2 -1 1 1 2 -1 1 8]
Candidate 2: Obj = 134.0474, schedule = [1 2 0 1 0 1 2 0 0 1 2 0 1 1 1 1 0 0 1 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 2 0 1 0 1 2 -1 2 1 1 0 0 1 1 1 1 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 2 1 1 -1 2 0 1 1 0 1 2 0 1 1 1 0 1 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 145
Iteration 19 completed in 5.49 seconds.
--- Iteration 20/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 1 0 1 0 1 2 0 0 1 1 2 0 1 0 1 2 -1 1 8]
Candidate 1: Obj = 132.0745, schedule = [1 1 2 0 0 1 2 0 0 1 1 1 2 0 0 1 1 1 0 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 1 2 -1 1 1 2 0 0 1 1 1 0 1 1 1 0 2 7]
Candidate 3: Obj = 123.1035, schedule = [2 1 2 0 0 1 1 1 0 1 2 1 0 0 1 1 1 0 1 8]
Candidate 4: Obj = 130.1952, schedule = [1 1 1 1 0 1 1 1 0 2 0 1 1 1 0 1 1 0 2 8]
Best objective value found so far: 114.9636
Total points evaluated: 150
Iteration 20 completed in 4.79 seconds.
--- Iteration 21/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 1 0 1 0 2 1 0 0 1 2 1 0 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 2 0 1 0 1 2 0 1 1 0 1 2 -1 1 1 2 0 0 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 2 0 2 0 0 1 2 1 0 1 1 1 0 1 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 2 1 1 0 0 1 1 1 1 1 0 2 0 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 1 2 0 0 2 1 0 1 1 0 1 1 0 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 155
Iteration 21 completed in 6.42 seconds.
--- Iteration 22/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 2 0 0 1 2 0 0 1 1 2 0 1 0 1 2 -1 2 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 1 1 0 1 2 1 -1 2 1 1 0 0 1 1 1 1 0 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 2 1 1 0 0 2 1 -1 2 1 0 2 0 1 0 1 0 2 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 2 1 0 1 1 1 0 2 0 1 1 0 1 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 1 1 1 0 1 1 0 0 1 2 0 2 -1 2 1 1 0 1 7]
Best objective value found so far: 114.9636
Total points evaluated: 160
Iteration 22 completed in 4.67 seconds.
--- Iteration 23/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 0 2 0 1 0 2 -1 2 0 2 1 0 1 0 2 0 1 8]
Candidate 1: Obj = 126.1819, schedule = [2 1 0 2 0 1 0 2 0 0 2 1 0 1 1 0 2 0 1 8]
Candidate 2: Obj = 131.9391, schedule = [2 1 0 2 0 1 0 2 0 1 0 2 1 0 1 0 2 0 0 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 1 0 2 -1 2 0 2 1 0 1 0 2 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 2 0 2 0 0 2 1 0 1 1 0 2 0 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 165
Iteration 23 completed in 4.24 seconds.
--- Iteration 24/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 0 1 1 0 1 2 -1 2 0 2 0 0 1 1 2 -1 2 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 2 0 1 1 0 2 0 0 2 0 2 1 0 0 1 2 -1 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 1 1 0 2 0 1 1 0 2 0 0 1 1 2 -1 2 8]
Candidate 3: Obj = 132.0435, schedule = [1 2 0 1 0 1 1 1 0 2 0 2 0 0 1 1 1 0 2 8]
Candidate 4: Obj = 128.0151, schedule = [1 2 0 1 0 1 2 0 0 1 1 2 0 0 1 1 1 1 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 170
Iteration 24 completed in 3.59 seconds.
--- Iteration 25/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 1 1 0 2 1 0 1 1 1 1 1 -1 1 2 1 -1 2 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 0 1 0 2 0 2 -1 2 1 0 2 0 0 2 1 -1 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 2 1 0 0 2 1 0 2 0 1 1 1 -1 2 8]
Candidate 3: Obj = 126.7964, schedule = [2 1 0 1 1 1 0 1 0 2 0 1 2 0 0 2 1 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 2 0 2 -1 2 0 1 2 0 0 2 1 0 1 8]
Best objective value found so far: 114.9636
Total points evaluated: 175
Iteration 25 completed in 4.78 seconds.
--- Optimization Finished ---
Total evaluations: 175
Best Objective Value Found: 114.96364292044709
Best U vector Found: [0 0 0 0 0 1 1 1 0 1 1 1 1 0 1 1 0 1 1 1]
--- Verification ---
Is the best U feasible? True
Resulting Y vector for best U: [2 1 1 1 1 1 1 0 1 1 1 1 0 1 1 0 2 0 1 7]
Objective value (recalculated): 114.96366.4.6.5 5. Experiment 2: CBO with Lower Confidence Bound (LCB) - Fixed Kappa
Applies the methodology from Deshwal et al. (2023) using LCB with fixed \(\kappa\).
# --- BO Helper Functions ---
# --- get_fitted_model function remains the same ---
def get_fitted_model(train_X_embedded_scaled, train_Y, m):
# ... (implementation is unchanged) ...
if train_Y.ndim > 1 and train_Y.shape[1] == 1: train_Y = train_Y.ravel()
kernel = ConstantKernel(1.0, constant_value_bounds=(1e-3, 1e3)) * \
Matern(length_scale=np.ones(m), length_scale_bounds=(1e-2, 1e2), nu=2.5) + \
WhiteKernel(noise_level=1e-10, # Small value for numerical stability
noise_level_bounds="fixed") # Bounds for noise optimization
gp_model = GaussianProcessRegressor(kernel=kernel, alpha=1e-10, n_restarts_optimizer=10, random_state=42)
gp_model.fit(train_X_embedded_scaled, train_Y)
return gp_model
def lower_confidence_bound(mu, sigma, kappa=2.576):
"""
Computes the Lower Confidence Bound (LCB) acquisition function.
Assumes maximization of this value guides the search (since mu is neg objective).
Higher LCB means lower predicted objective or lower penalty for uncertainty.
mu, sigma: Predicted mean and standard deviation (NumPy arrays).
kappa: Controls the balance between exploitation (high mu -> low original objective)
and exploration (low sigma).
"""
# Ensure sigma is non-negative
sigma = np.maximum(sigma, 0)
return mu - kappa * sigma # <<< Sign flipped from UCB
def optimize_acqf_discrete_via_embedding(gp_model, scaler, dictionary_A, T, q, num_candidates, kappa):
"""
Optimizes LCB acquisition function by sampling random binary candidates,
embedding, SCALING, predicting with GP, and calculating LCB.
Selects the top q candidates based on LCB.
Returns candidates as a numpy array (q x T).
"""
m = dictionary_A.shape[0]
# 1. Generate Random Binary Candidates
candidate_u_vectors_np = np.random.randint(0, 2, size=(num_candidates, T))
# 2. Embed the Candidates
embedded_candidates_np = embed_batch(candidate_u_vectors_np, dictionary_A)
# 3. Scale the Embedded Candidates
embedded_candidates_scaled_np = scaler.transform(embedded_candidates_np)
# 4. Predict Mean and Std Dev using the GP Model
mu, std = gp_model.predict(embedded_candidates_scaled_np, return_std=True)
# 5. Calculate Acquisition Function (Lower Confidence Bound) <<< CHANGED HERE
acq_values = lower_confidence_bound(mu, std, kappa=kappa) # Use LCB
# 6. Select Top Candidates (based on highest LCB) <<< COMMENT UPDATED
# We maximize LCB = mu - kappa*sigma, where mu is neg_objective
top_indices = np.argsort(acq_values)[-q:]
top_indices = top_indices[::-1] # Ensure descending order of LCB
return candidate_u_vectors_np[top_indices, :]
# --- BO Loop ---
# Parameters
KAPPA = 2.576 # Exploration parameter for LCB. Adjust as needed.
N_INITIAL = 49 # The initial schedule will always be included
N_ITERATIONS = 25
BATCH_SIZE_q = 5
NUM_CANDIDATES_Acqf = T*5*1024 # Might need more for higher T
m = math.ceil(T/4) # Dimension of the embedding space
# Store evaluated points (using NumPy arrays)
evaluated_U_np_list = [] # List to store evaluated U vectors (binary)
evaluated_f_vals = [] # List to store raw objective values (lower is better)
train_Y_list = [] # List to store NEGATED objective values for GP (higher is better)
# 1. Initialization
for U_init in initial_candidates:
f_val = evaluate_objective(U_init, X, v_star, convolutions, d, w)
neg_f_val = -f_val
evaluated_U_np_list.append(U_init)
evaluated_f_vals.append(f_val)
train_Y_list.append(neg_f_val)
# Convert lists to NumPy arrays for GP fitting
train_Y = np.array(train_Y_list).reshape(-1, 1) # Keep as column vector initially
best_obj_so_far = min(evaluated_f_vals) if evaluated_f_vals else float('inf')
initial_best_obj_so_far_lcb = best_obj_so_far # Saved for reporting results
print(f"Initial best objective value: {best_obj_so_far}")
if not np.isfinite(best_obj_so_far):
print("Warning: Initial best objective is infinite, possibly all initial points were infeasible.")
# 2. BO Iterations
for iteration in range(N_ITERATIONS):
start_time = time.time()
print(f"\n--- Iteration {iteration + 1}/{N_ITERATIONS} ---")
# a. Generate dictionary A (remains the same)
current_dictionary_A = generate_diverse_random_dictionary(T, m)
# b. Embed ALL evaluated U vectors (remains the same)
if not evaluated_U_np_list: continue
evaluated_U_np_array = np.array(evaluated_U_np_list)
embedded_train_X = embed_batch(evaluated_U_np_array, current_dictionary_A)
# c. Scale the embedded training data (remains the same)
scaler = MinMaxScaler()
if embedded_train_X.shape[0] > 0: embedded_train_X_scaled = scaler.fit_transform(embedded_train_X)
else: embedded_train_X_scaled = embedded_train_X
# Ensure train_Y is NumPy array
train_Y_for_fit = np.array(train_Y_list)
# d. Fit GP Model (remains the same)
print("Fitting GP model...")
if embedded_train_X_scaled.shape[0] > 0 and train_Y_for_fit.shape[0] == embedded_train_X_scaled.shape[0]:
gp_model = get_fitted_model(embedded_train_X_scaled, train_Y_for_fit, m)
print("GP model fitted.")
else:
print("Warning: Not enough data or data mismatch to fit GP model. Skipping iteration.")
continue
# e. Determine current best neg value (useful for tracking, not directly used in LCB calculation)
current_best_neg_f_val = np.max(train_Y_for_fit) if train_Y_for_fit.size > 0 else -float('inf')
if current_best_neg_f_val <= -LARGE_PENALTY / 2 and np.isfinite(current_best_neg_f_val):
print(f"Warning: Current best NEGATIVE objective value ({current_best_neg_f_val:.2f}) is very low (likely from penalties).")
# f. Optimize Acquisition Function (LCB) <<< MODIFIED CALL & COMMENT
print("Optimizing acquisition function (LCB)...") # Comment updated
next_U_candidates_np = optimize_acqf_discrete_via_embedding(
gp_model=gp_model,
scaler=scaler,
dictionary_A=current_dictionary_A,
T=T,
q=BATCH_SIZE_q,
num_candidates=NUM_CANDIDATES_Acqf,
kappa=KAPPA # Pass kappa
)
print(f"Selected {next_U_candidates_np.shape[0]} candidate(s).")
# g. Evaluate Objective (remains the same)
newly_evaluated_U = []; newly_evaluated_f = []; newly_evaluated_neg_f = []
for i in range(next_U_candidates_np.shape[0]):
next_U = next_U_candidates_np[i, :]
already_evaluated = any(np.array_equal(next_U, u) for u in evaluated_U_np_list)
if already_evaluated: print(f" Candidate {i} was already evaluated. Skipping."); continue
next_f = evaluate_objective(next_U, X, v_star, convolutions, d, w)
next_neg_f = -next_f
V_sum = np.sum(v_star[next_U == 1, :], axis=0)
Y = X + V_sum
print(f" Candidate {i}: Obj = {next_f:.4f}, schedule = {Y}")
newly_evaluated_U.append(next_U); newly_evaluated_f.append(next_f); newly_evaluated_neg_f.append(next_neg_f)
if next_f < best_obj_so_far: best_obj_so_far = next_f
# h. Augment Dataset (remains the same)
evaluated_U_np_list.extend(newly_evaluated_U); evaluated_f_vals.extend(newly_evaluated_f); train_Y_list.extend(newly_evaluated_neg_f)
train_Y = np.array(train_Y_list).reshape(-1, 1)
iter_time = time.time() - start_time
print(f"Best objective value found so far: {best_obj_so_far:.4f}")
print(f"Total points evaluated: {len(evaluated_f_vals)}")
print(f"Iteration {iteration + 1} completed in {iter_time:.2f} seconds.")
# --- Results ---
print("\n--- Optimization Finished ---")
if not evaluated_f_vals: print("No points were successfully evaluated.")
else:
final_best_idx_lcb = np.argmin(evaluated_f_vals)
final_best_U_lcb = evaluated_U_np_list[final_best_idx_lcb]
final_best_f_lcb = evaluated_f_vals[final_best_idx_lcb]
nr_evaluated_f_vals_lcb = len(evaluated_f_vals) # Saved for reporting results
print(f"Total evaluations: {nr_evaluated_f_vals_lcb}")
print(f"Best Objective Value Found: {final_best_f_lcb}")
print(f"Best U vector Found: {final_best_U_lcb}")
# Verification
V_sum_best = np.sum(v_star[final_best_U_lcb == 1, :], axis=0)
Y_best_lcb = X + V_sum_best
is_feasible = np.all(Y_best_lcb >= 0)
recalculated_obj = LARGE_PENALTY
if is_feasible:
ewt, esp = calculate_objective_serv_time_lookup(Y_best_lcb, d, convolutions)
recalculated_obj = w * ewt + (1 - w) * esp
print(f"\n--- Verification ---")
print(f"Is the best U feasible? {is_feasible}")
if is_feasible:
print(f"Resulting Y vector for best U: {Y_best_lcb}")
print(f"Objective value (recalculated): {recalculated_obj:.4f}")
if not np.isclose(final_best_f_lcb, recalculated_obj): print(f"Warning: Stored best objective ({final_best_f:.4f}) does not match recalculation ({recalculated_obj:.4f})!")
elif final_best_f_lcb < LARGE_PENALTY: print(f"Warning: Best objective ({final_best_f_lcb:.4f}) is not the penalty value, but feasibility check failed."); print(f"Resulting Y vector (infeasible): {Y_best_lcb}")
else: print("Best solution found corresponds to an infeasible penalty value.")Initial best objective value: 118.87100166955729
--- Iteration 1/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 1 1 1 1 1 0 0 2 1 1 0 1 0 1 2 -1 1 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 1 0 2 -1 1 1 2 0 0 1 1 1 1 1 8]
Candidate 2: Obj = 133.2394, schedule = [1 2 0 2 0 1 1 0 0 2 1 0 1 0 2 0 1 1 0 9]
Candidate 3: Obj = 134.0944, schedule = [2 1 1 0 0 2 0 1 0 2 0 2 0 0 2 1 0 1 0 9]
Candidate 4: Obj = 135.2522, schedule = [2 1 1 1 0 1 1 0 0 1 2 1 0 0 1 1 1 0 1 9]
Best objective value found so far: 118.8710
Total points evaluated: 55
Iteration 1 completed in 4.63 seconds.
--- Iteration 2/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 129.3169, schedule = [1 2 1 1 0 0 1 1 0 2 1 1 0 0 1 2 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 1 0 2 0 0 1 1 0 2 1 1 1 -1 1 2 1 0 1 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 2 0 2 -1 2 1 1 1 0 0 1 2 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 2 0 2 0 0 1 2 -1 2 1 1 1 -1 1 2 1 0 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 1 0 2 0 0 1 2 -1 2 1 1 1 -1 1 1 2 0 1 7]
Best objective value found so far: 118.8710
Total points evaluated: 60
Iteration 2 completed in 4.04 seconds.
--- Iteration 3/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 2 1 0 0 2 1 -1 2 0 1 1 1 0 2 0 1 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 0 2 -1 2 1 1 1 0 1 1 0 1 0 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 0 2 1 -1 2 1 0 1 1 1 1 0 1 1 1 0 1 0 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 1 1 0 1 1 1 0 2 0 1 2 -1 2 0 2 0 1 7]
Candidate 4: Obj = 117.6316, schedule = [2 1 1 1 0 1 2 0 0 1 2 0 1 0 1 2 1 0 1 7]
Best objective value found so far: 117.6316
Total points evaluated: 65
Iteration 3 completed in 4.27 seconds.
--- Iteration 4/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 1 0 1 0 2 1 1 0 1 0 2 0 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 1 0 1 0 2 1 0 1 1 1 1 0 0 8]
Candidate 2: Obj = 117.4455, schedule = [2 1 1 2 0 1 1 0 1 0 2 0 2 0 1 0 2 0 1 7]
Candidate 3: Obj = 118.6866, schedule = [3 0 1 2 0 1 1 0 1 0 2 1 0 1 0 2 0 1 1 7]
Candidate 4: Obj = 115.6020, schedule = [2 1 2 0 1 1 0 1 1 0 2 0 2 0 1 1 1 0 1 7]
Best objective value found so far: 115.6020
Total points evaluated: 70
Iteration 4 completed in 4.70 seconds.
--- Iteration 5/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 0 2 1 0 0 2 1 -1 1 1 2 1 0 1 0 2 0 1 7]
Candidate 1: Obj = 120.3275, schedule = [3 1 1 1 0 1 1 1 0 0 1 1 1 0 2 1 0 1 1 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 1 1 0 1 1 1 1 -1 2 0 2 1 -1 2 1 0 1 0 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 0 2 1 -1 2 1 1 1 -1 2 0 1 0 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 1 2 1 0 1 1 0 0 2 0 2 1 -1 2 0 2 0 1 8]
Best objective value found so far: 115.6020
Total points evaluated: 75
Iteration 5 completed in 3.88 seconds.
--- Iteration 6/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 1 1 -1 2 1 0 1 0 1 1 1 0 1 2 1 0 1 8]
Candidate 1: Obj = 128.5334, schedule = [2 1 1 0 0 1 2 1 0 0 1 1 1 0 1 2 1 0 1 8]
Candidate 2: Obj = 127.5458, schedule = [2 0 2 0 0 2 1 0 1 1 0 1 1 0 1 2 1 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 0 2 1 -1 1 2 0 1 1 1 0 2 -1 1 2 1 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 1 1 -1 2 1 0 1 1 0 1 1 0 1 2 1 -1 1 9]
Best objective value found so far: 115.6020
Total points evaluated: 80
Iteration 6 completed in 3.95 seconds.
--- Iteration 7/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 1 0 1 1 1 0 1 1 0 1 2 -1 2 0 1 1 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 1 2 0 1 1 0 2 -1 2 0 1 2 0 0 2 1 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 0 2 0 0 2 1 1 -1 2 0 1 2 -1 1 2 1 0 1 8]
Candidate 3: Obj = 125.7593, schedule = [2 0 2 0 0 2 0 1 1 1 0 1 2 0 1 0 2 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 2 0 2 0 1 0 1 2 -1 1 2 1 0 1 8]
Best objective value found so far: 115.6020
Total points evaluated: 85
Iteration 7 completed in 4.11 seconds.
--- Iteration 8/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 123.2298, schedule = [2 1 2 0 1 0 2 0 1 0 1 1 2 0 0 1 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 2 0 1 0 2 0 1 0 2 0 2 -1 1 1 1 0 1 8]
Candidate 2: Obj = 123.1074, schedule = [2 1 2 0 0 1 2 0 1 1 0 1 2 0 0 1 1 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 2 1 0 0 1 2 0 1 0 1 2 1 -1 1 1 1 0 1 8]
Candidate 4: Obj = 124.3234, schedule = [2 1 2 0 0 2 1 0 1 0 2 1 0 0 1 1 1 0 1 8]
Best objective value found so far: 115.6020
Total points evaluated: 90
Iteration 8 completed in 3.96 seconds.
--- Iteration 9/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 130.8021, schedule = [1 2 0 2 0 1 1 1 0 1 1 1 1 0 0 1 1 1 0 9]
Candidate 1: Obj = 130.5450, schedule = [2 1 0 2 0 1 1 1 0 1 0 1 2 0 1 0 1 1 0 9]
Candidate 2: Obj = 124.2772, schedule = [2 1 0 2 0 1 1 1 0 1 0 1 1 1 1 0 1 1 1 8]
Candidate 3: Obj = 126.1602, schedule = [1 2 0 2 0 1 1 1 0 1 0 2 0 1 0 1 2 0 1 8]
Candidate 4: Obj = 132.6468, schedule = [2 1 0 2 0 1 1 1 0 0 2 1 1 0 0 1 1 1 0 9]
Best objective value found so far: 115.6020
Total points evaluated: 95
Iteration 9 completed in 4.18 seconds.
--- Iteration 10/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 0 2 -1 2 1 0 1 1 0 2 1 0 1 1 0 1 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 2 0 2 0 0 1 2 1 0 1 0 2 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 2 0 2 0 0 1 2 1 0 1 0 2 0 1 8]
Candidate 3: Obj = 128.4867, schedule = [1 1 2 0 0 2 1 1 0 0 2 1 1 0 0 2 1 0 1 8]
Candidate 4: Obj = 130.2438, schedule = [1 1 1 1 0 2 1 1 0 1 1 1 1 0 1 1 1 0 0 9]
Best objective value found so far: 115.6020
Total points evaluated: 100
Iteration 10 completed in 5.25 seconds.
--- Iteration 11/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 130.5600, schedule = [2 0 1 1 1 0 2 0 1 1 0 1 1 1 0 2 0 0 2 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 1 2 0 1 1 1 1 0 0 1 1 1 0 2 8]
Candidate 2: Obj = 134.2707, schedule = [2 0 1 1 0 2 1 0 1 1 1 1 0 0 1 2 0 0 1 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 0 1 0 1 2 0 1 0 1 2 1 -1 1 2 1 -1 1 9]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 1 2 0 1 0 1 2 1 0 0 2 0 0 2 8]
Best objective value found so far: 115.6020
Total points evaluated: 105
Iteration 11 completed in 4.73 seconds.
--- Iteration 12/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 2 1 1 0 1 1 1 1 0 1 1 1 0 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 0 2 1 -1 2 1 0 2 0 1 1 1 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 1 0 2 0 1 1 1 0 1 0 2 1 -1 2 1 1 0 1 7]
Candidate 3: Obj = 125.1870, schedule = [1 2 0 2 0 0 2 1 0 1 1 1 1 0 1 0 2 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 2 0 1 1 1 -1 2 1 0 2 0 1 1 1 0 1 8]
Best objective value found so far: 115.6020
Total points evaluated: 110
Iteration 12 completed in 5.29 seconds.
--- Iteration 13/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 131.4070, schedule = [2 1 1 0 1 0 1 1 0 1 1 2 0 1 0 1 1 1 0 9]
Candidate 1: Obj = 122.3589, schedule = [2 2 1 0 1 0 1 2 0 0 1 2 0 1 0 1 1 1 0 8]
Candidate 2: Obj = 131.4070, schedule = [2 1 1 0 1 0 1 1 0 1 1 2 0 1 0 1 1 1 0 9]
Candidate 3: Obj = 131.4662, schedule = [2 1 1 0 1 0 1 1 1 0 1 2 0 0 1 1 1 1 0 9]
Candidate 4: Obj = 122.4985, schedule = [3 1 0 1 1 0 1 2 0 0 1 2 0 1 0 1 1 0 2 7]
Best objective value found so far: 115.6020
Total points evaluated: 115
Iteration 13 completed in 8.04 seconds.
--- Iteration 14/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 122.3147, schedule = [3 0 2 0 0 1 2 0 0 1 2 0 1 1 1 1 0 1 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 2 1 -1 1 1 2 -1 1 1 1 1 0 2 1 0 1 0 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 2 0 0 1 2 1 -1 2 1 0 1 0 2 0 1 1 0 8]
Candidate 3: Obj = 115.1129, schedule = [2 1 1 1 0 1 2 0 0 1 2 0 1 1 1 1 1 0 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 2 1 -1 1 2 0 0 1 1 2 0 0 2 0 1 1 0 8]
Best objective value found so far: 115.1129
Total points evaluated: 120
Iteration 14 completed in 7.87 seconds.
--- Iteration 15/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 123.3545, schedule = [2 1 2 0 0 2 1 1 0 1 0 1 1 0 2 0 2 0 0 8]
Candidate 1: Obj = 118.9808, schedule = [2 1 2 0 0 2 1 1 0 0 1 1 1 0 2 0 2 0 1 7]
Candidate 2: Obj = 124.7961, schedule = [3 1 0 2 0 0 2 0 0 2 0 2 0 1 0 2 0 0 2 7]
Candidate 3: Obj = 122.3054, schedule = [2 2 0 2 0 0 1 1 1 1 1 1 1 0 0 1 1 0 2 7]
Candidate 4: Obj = 133.7135, schedule = [2 0 2 0 0 2 1 0 0 2 0 1 1 0 2 1 0 1 0 9]
Best objective value found so far: 115.1129
Total points evaluated: 125
Iteration 15 completed in 6.16 seconds.
--- Iteration 16/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 124.3527, schedule = [3 1 0 2 0 0 2 0 1 1 0 2 0 0 1 1 2 0 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 1 2 -1 1 2 0 0 1 2 1 1 -1 2 0 2 -1 2 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 2 0 1 1 0 1 1 1 0 2 1 0 0 2 0 2 -1 1 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 0 2 -1 2 1 1 0 0 1 1 2 0 0 8]
Candidate 4: Obj = 132.6420, schedule = [2 0 1 2 0 0 1 1 0 1 2 0 1 0 1 2 0 0 2 8]
Best objective value found so far: 115.1129
Total points evaluated: 130
Iteration 16 completed in 6.82 seconds.
--- Iteration 17/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 2 1 0 1 1 0 1 0 1 2 0 1 0 1 1 2 -1 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 1 0 1 1 1 0 1 0 1 2 0 1 0 1 1 2 -1 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 2 0 1 1 0 1 1 1 0 2 1 1 0 0 1 2 -1 1 8]
Candidate 3: Obj = 123.5410, schedule = [2 2 0 1 1 1 1 0 0 1 2 0 1 1 0 1 1 0 1 8]
Candidate 4: Obj = 122.4941, schedule = [2 2 0 1 1 0 1 1 1 1 0 1 2 0 0 1 1 0 1 8]
Best objective value found so far: 115.1129
Total points evaluated: 135
Iteration 17 completed in 7.43 seconds.
--- Iteration 18/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 2 1 0 0 1 2 1 0 1 0 2 1 -1 2 0 1 1 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 2 0 2 -1 2 1 1 0 1 0 2 0 1 0 1 1 1 0 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 2 1 1 -1 2 1 1 0 0 1 1 1 1 0 2 0 1 0 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 2 1 1 -1 2 0 2 0 1 0 1 1 0 1 2 0 1 0 9]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 1 0 2 -1 2 1 1 0 0 1 1 1 0 2 1 0 1 0 8]
Best objective value found so far: 115.1129
Total points evaluated: 140
Iteration 18 completed in 8.02 seconds.
--- Iteration 19/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 1 -1 2 1 1 0 0 2 1 1 -1 2 1 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 2 0 1 0 1 1 2 -1 1 1 1 1 0 1 1 1 1 1 7]
Candidate 2: Obj = 132.4621, schedule = [2 1 1 0 1 1 0 2 0 0 2 1 1 0 1 1 0 0 2 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 1 1 2 -1 2 1 0 1 0 2 0 1 0 1 8]
Candidate 4: Obj = 133.4902, schedule = [1 1 1 1 0 1 2 0 0 2 0 1 1 0 2 0 1 0 1 9]
Best objective value found so far: 115.1129
Total points evaluated: 145
Iteration 19 completed in 6.84 seconds.
--- Iteration 20/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 2 0 1 0 1 0 2 1 1 1 -1 2 1 0 1 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 1 1 2 0 1 0 2 -1 2 0 2 0 1 0 1 1 0 2 8]
Candidate 2: Obj = 126.6556, schedule = [1 2 0 1 0 1 1 2 0 1 1 0 1 0 2 1 0 1 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 0 0 2 1 1 0 0 1 1 1 1 0 1 2 -1 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 2 1 0 0 1 2 -1 2 0 1 1 0 2 1 0 1 1 7]
Best objective value found so far: 115.1129
Total points evaluated: 150
Iteration 20 completed in 8.16 seconds.
--- Iteration 21/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 123.3261, schedule = [2 1 2 1 0 1 0 2 0 0 1 1 2 0 1 0 2 0 0 8]
Candidate 1: Obj = 123.2211, schedule = [2 1 2 1 0 1 0 2 0 0 1 2 0 1 1 0 2 0 0 8]
Candidate 2: Obj = 122.8762, schedule = [2 1 2 1 0 1 1 1 0 0 1 1 2 0 1 0 1 1 0 8]
Candidate 3: Obj = 132.0113, schedule = [1 1 2 1 0 1 0 2 0 0 1 1 2 0 1 0 1 0 2 8]
Candidate 4: Obj = 118.8922, schedule = [2 2 0 2 0 1 1 0 0 1 1 2 1 0 1 0 1 1 1 7]
Best objective value found so far: 115.1129
Total points evaluated: 155
Iteration 21 completed in 7.71 seconds.
--- Iteration 22/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 124.9031, schedule = [3 0 2 1 0 0 2 1 0 1 1 0 1 0 2 0 2 0 0 8]
Candidate 1: Obj = 120.0952, schedule = [3 0 2 1 0 0 2 1 0 1 0 1 1 0 2 0 2 0 1 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 2 1 0 0 2 0 1 1 1 0 2 -1 1 2 1 0 0 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 1 2 0 1 1 1 0 1 1 0 2 -1 2 0 1 1 0 8]
Candidate 4: Obj = 119.4426, schedule = [2 1 2 1 0 0 2 1 0 1 1 0 1 0 2 0 2 0 1 7]
Best objective value found so far: 115.1129
Total points evaluated: 160
Iteration 22 completed in 6.18 seconds.
--- Iteration 23/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 0 1 1 1 0 1 0 2 0 2 -1 2 0 2 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 2 0 2 -1 2 1 1 0 1 1 1 0 1 1 0 2 0 1 7]
Candidate 2: Obj = 123.1513, schedule = [1 1 1 1 1 1 0 2 0 1 0 2 0 1 1 1 0 1 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 0 1 1 1 1 0 1 1 1 1 1 -1 2 0 2 -1 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 0 1 0 2 1 1 0 1 1 0 1 1 1 1 1 -1 2 8]
Best objective value found so far: 115.1129
Total points evaluated: 165
Iteration 23 completed in 4.94 seconds.
--- Iteration 24/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 125.8798, schedule = [2 1 1 1 0 0 1 1 1 1 1 1 1 0 0 2 1 0 1 8]
Candidate 1: Obj = 127.9569, schedule = [2 1 1 1 0 0 2 0 1 0 2 1 1 0 0 2 1 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 0 2 0 1 1 1 1 1 0 0 2 1 -1 2 8]
Candidate 3: Obj = 127.9569, schedule = [2 1 1 1 0 0 2 0 1 0 2 1 1 0 0 2 1 0 1 8]
Candidate 4: Obj = 126.5368, schedule = [2 1 1 1 0 0 2 0 1 1 1 1 1 0 0 1 1 1 1 8]
Best objective value found so far: 115.1129
Total points evaluated: 170
Iteration 24 completed in 6.99 seconds.
--- Iteration 25/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 2 1 0 0 2 1 -1 2 0 2 1 -1 1 2 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 2 0 0 2 1 1 0 1 1 1 1 -1 1 1 2 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 1 1 0 1 1 1 1 -1 1 2 1 0 1 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 1 2 0 1 0 2 1 0 1 1 1 1 -1 1 2 1 0 1 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 1 2 1 -1 2 1 1 0 1 1 1 1 -1 1 2 1 0 0 9]
Best objective value found so far: 115.1129
Total points evaluated: 175
Iteration 25 completed in 6.32 seconds.
--- Optimization Finished ---
Total evaluations: 175
Best Objective Value Found: 115.11287490459011
Best U vector Found: [0 0 0 0 0 0 0 1 0 0 0 1 0 0 1 1 1 1 1 1]
--- Verification ---
Is the best U feasible? True
Resulting Y vector for best U: [2 1 1 1 0 1 2 0 0 1 2 0 1 1 1 1 1 0 1 7]
Objective value (recalculated): 115.11296.4.6.6 6. Experiment 3: CBO with LCB and wavelet dictionary.
# --- BO Helper Functions ---
# --- get_fitted_model function remains the same ---
def get_fitted_model(train_X_embedded_scaled, train_Y, m):
# ... (implementation is unchanged) ...
if train_Y.ndim > 1 and train_Y.shape[1] == 1: train_Y = train_Y.ravel()
kernel = ConstantKernel(1.0, constant_value_bounds=(1e-3, 1e3)) * \
Matern(length_scale=np.ones(m), length_scale_bounds=(1e-2, 1e2), nu=2.5) + \
WhiteKernel(noise_level=1e-10, # Small value for numerical stability
noise_level_bounds="fixed") # Bounds for noise optimization
gp_model = GaussianProcessRegressor(kernel=kernel, alpha=1e-10, n_restarts_optimizer=10, random_state=42)
gp_model.fit(train_X_embedded_scaled, train_Y)
return gp_model
def lower_confidence_bound(mu, sigma, kappa=2.576):
"""
Computes the Lower Confidence Bound (LCB) acquisition function.
Assumes maximization of this value guides the search (since mu is neg objective).
Higher LCB means lower predicted objective or lower penalty for uncertainty.
mu, sigma: Predicted mean and standard deviation (NumPy arrays).
kappa: Controls the balance between exploitation (high mu -> low original objective)
and exploration (low sigma).
"""
# Ensure sigma is non-negative
sigma = np.maximum(sigma, 0)
return mu - kappa * sigma # <<< Sign flipped from UCB
def optimize_acqf_discrete_via_embedding(gp_model, scaler, dictionary_A, T, q, num_candidates, kappa):
"""
Optimizes LCB acquisition function by sampling random binary candidates,
embedding, SCALING, predicting with GP, and calculating LCB.
Selects the top q candidates based on LCB.
Returns candidates as a numpy array (q x T).
"""
m = dictionary_A.shape[0]
# 1. Generate Random Binary Candidates
candidate_u_vectors_np = np.random.randint(0, 2, size=(num_candidates, T))
# 2. Embed the Candidates
embedded_candidates_np = embed_batch(candidate_u_vectors_np, dictionary_A)
# 3. Scale the Embedded Candidates
embedded_candidates_scaled_np = scaler.transform(embedded_candidates_np)
# 4. Predict Mean and Std Dev using the GP Model
mu, std = gp_model.predict(embedded_candidates_scaled_np, return_std=True)
# 5. Calculate Acquisition Function (Lower Confidence Bound) <<< CHANGED HERE
acq_values = lower_confidence_bound(mu, std, kappa=kappa) # Use LCB
# 6. Select Top Candidates (based on highest LCB) <<< COMMENT UPDATED
# We maximize LCB = mu - kappa*sigma, where mu is neg_objective
top_indices = np.argsort(acq_values)[-q:]
top_indices = top_indices[::-1] # Ensure descending order of LCB
return candidate_u_vectors_np[top_indices, :]
# --- BO Loop ---
# Parameters
KAPPA = 2.576 # Exploration parameter for LCB. Adjust as needed.
N_INITIAL = 49 # The initial schedule will always be included
N_ITERATIONS = 25
BATCH_SIZE_q = 5
NUM_CANDIDATES_Acqf = T*5*1024 # Might need more for higher T
m = math.ceil(T/4) # Dimension of the embedding space
# Store evaluated points (using NumPy arrays)
evaluated_U_np_list = [] # List to store evaluated U vectors (binary)
evaluated_f_vals = [] # List to store raw objective values (lower is better)
train_Y_list = [] # List to store NEGATED objective values for GP (higher is better)
# 1. Initialization
for U_init in initial_candidates:
f_val = evaluate_objective(U_init, X, v_star, convolutions, d, w)
neg_f_val = -f_val
evaluated_U_np_list.append(U_init)
evaluated_f_vals.append(f_val)
train_Y_list.append(neg_f_val)
# Convert lists to NumPy arrays for GP fitting
train_Y = np.array(train_Y_list).reshape(-1, 1) # Keep as column vector initially
best_obj_so_far = min(evaluated_f_vals) if evaluated_f_vals else float('inf')
initial_best_obj_so_far_lcb = best_obj_so_far # Saved for reporting results
print(f"Initial best objective value: {best_obj_so_far}")
if not np.isfinite(best_obj_so_far):
print("Warning: Initial best objective is infinite, possibly all initial points were infeasible.")
# 2. BO Iterations
for iteration in range(N_ITERATIONS):
start_time = time.time()
print(f"\n--- Iteration {iteration + 1}/{N_ITERATIONS} ---")
# a. Generate dictionary A (remains the same)
current_dictionary_A = generate_wavelet_dictionary(T, m)
# b. Embed ALL evaluated U vectors (remains the same)
if not evaluated_U_np_list: continue
evaluated_U_np_array = np.array(evaluated_U_np_list)
embedded_train_X = embed_batch(evaluated_U_np_array, current_dictionary_A)
# c. Scale the embedded training data (remains the same)
scaler = MinMaxScaler()
if embedded_train_X.shape[0] > 0: embedded_train_X_scaled = scaler.fit_transform(embedded_train_X)
else: embedded_train_X_scaled = embedded_train_X
# Ensure train_Y is NumPy array
train_Y_for_fit = np.array(train_Y_list)
# d. Fit GP Model (remains the same)
print("Fitting GP model...")
if embedded_train_X_scaled.shape[0] > 0 and train_Y_for_fit.shape[0] == embedded_train_X_scaled.shape[0]:
gp_model = get_fitted_model(embedded_train_X_scaled, train_Y_for_fit, m)
print("GP model fitted.")
else:
print("Warning: Not enough data or data mismatch to fit GP model. Skipping iteration.")
continue
# e. Determine current best neg value (useful for tracking, not directly used in LCB calculation)
current_best_neg_f_val = np.max(train_Y_for_fit) if train_Y_for_fit.size > 0 else -float('inf')
if current_best_neg_f_val <= -LARGE_PENALTY / 2 and np.isfinite(current_best_neg_f_val):
print(f"Warning: Current best NEGATIVE objective value ({current_best_neg_f_val:.2f}) is very low (likely from penalties).")
# f. Optimize Acquisition Function (LCB) <<< MODIFIED CALL & COMMENT
print("Optimizing acquisition function (LCB)...") # Comment updated
next_U_candidates_np = optimize_acqf_discrete_via_embedding(
gp_model=gp_model,
scaler=scaler,
dictionary_A=current_dictionary_A,
T=T,
q=BATCH_SIZE_q,
num_candidates=NUM_CANDIDATES_Acqf,
kappa=KAPPA # Pass kappa
)
print(f"Selected {next_U_candidates_np.shape[0]} candidate(s).")
# g. Evaluate Objective (remains the same)
newly_evaluated_U = []; newly_evaluated_f = []; newly_evaluated_neg_f = []
for i in range(next_U_candidates_np.shape[0]):
next_U = next_U_candidates_np[i, :]
already_evaluated = any(np.array_equal(next_U, u) for u in evaluated_U_np_list)
if already_evaluated: print(f" Candidate {i} was already evaluated. Skipping."); continue
next_f = evaluate_objective(next_U, X, v_star, convolutions, d, w)
next_neg_f = -next_f
V_sum = np.sum(v_star[next_U == 1, :], axis=0)
Y = X + V_sum
print(f" Candidate {i}: Obj = {next_f:.4f}, schedule = {Y}")
newly_evaluated_U.append(next_U); newly_evaluated_f.append(next_f); newly_evaluated_neg_f.append(next_neg_f)
if next_f < best_obj_so_far: best_obj_so_far = next_f
# h. Augment Dataset (remains the same)
evaluated_U_np_list.extend(newly_evaluated_U); evaluated_f_vals.extend(newly_evaluated_f); train_Y_list.extend(newly_evaluated_neg_f)
train_Y = np.array(train_Y_list).reshape(-1, 1)
iter_time = time.time() - start_time
print(f"Best objective value found so far: {best_obj_so_far:.4f}")
print(f"Total points evaluated: {len(evaluated_f_vals)}")
print(f"Iteration {iteration + 1} completed in {iter_time:.2f} seconds.")
# --- Results ---
print("\n--- Optimization Finished ---")
if not evaluated_f_vals: print("No points were successfully evaluated.")
else:
final_best_idx_lcb = np.argmin(evaluated_f_vals)
final_best_U_lcb = evaluated_U_np_list[final_best_idx_lcb]
final_best_f_lcb = evaluated_f_vals[final_best_idx_lcb]
nr_evaluated_f_vals_lcb = len(evaluated_f_vals) # Saved for reporting results
print(f"Total evaluations: {nr_evaluated_f_vals_lcb}")
print(f"Best Objective Value Found: {final_best_f_lcb}")
print(f"Best U vector Found: {final_best_U_lcb}")
# Verification
V_sum_best = np.sum(v_star[final_best_U_lcb == 1, :], axis=0)
Y_best_lcb = X + V_sum_best
is_feasible = np.all(Y_best_lcb >= 0)
recalculated_obj = LARGE_PENALTY
if is_feasible:
ewt, esp = calculate_objective_serv_time_lookup(Y_best_lcb, d, convolutions)
recalculated_obj = w * ewt + (1 - w) * esp
print(f"\n--- Verification ---")
print(f"Is the best U feasible? {is_feasible}")
if is_feasible:
print(f"Resulting Y vector for best U: {Y_best_lcb}")
print(f"Objective value (recalculated): {recalculated_obj:.4f}")
if not np.isclose(final_best_f_lcb, recalculated_obj): print(f"Warning: Stored best objective ({final_best_f:.4f}) does not match recalculation ({recalculated_obj:.4f})!")
elif final_best_f_lcb < LARGE_PENALTY: print(f"Warning: Best objective ({final_best_f_lcb:.4f}) is not the penalty value, but feasibility check failed."); print(f"Resulting Y vector (infeasible): {Y_best_lcb}")
else: print("Best solution found corresponds to an infeasible penalty value.")Initial best objective value: 118.87100166955729
--- Iteration 1/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 2 1 0 0 2 0 2 1 0 0 1 1 1 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 1 1 -1 2 0 1 0 1 1 2 1 0 1 1 0 0 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 2 0 1 1 1 1 1 1 -1 2 0 1 1 1 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 2 0 0 1 1 1 0 2 0 1 1 1 0 2 1 -1 1 8]
Candidate 4: Obj = 123.5187, schedule = [2 1 1 1 0 2 1 0 0 1 1 2 0 1 1 1 0 0 1 8]
Best objective value found so far: 118.8710
Total points evaluated: 55
Iteration 1 completed in 3.57 seconds.
--- Iteration 2/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 130.8475, schedule = [1 2 0 1 1 1 0 2 0 1 1 0 1 1 0 2 1 0 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 2 0 1 1 1 0 2 -1 2 0 2 0 1 0 2 1 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 2 1 0 1 1 1 1 0 1 0 2 0 0 1 2 1 -1 1 9]
Candidate 3: Obj = 121.4522, schedule = [2 1 2 0 1 1 1 0 1 1 0 2 0 1 1 1 1 0 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 1 2 0 0 2 1 1 0 1 0 2 1 -1 1 2 0 1 1 8]
Best objective value found so far: 118.8710
Total points evaluated: 60
Iteration 2 completed in 3.65 seconds.
--- Iteration 3/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 131.7503, schedule = [1 2 0 2 0 1 0 2 0 1 1 1 0 1 0 2 1 0 0 9]
Candidate 1: Obj = 124.4007, schedule = [3 1 0 1 1 1 1 1 0 1 0 2 1 0 0 2 1 0 0 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 2 0 2 0 1 1 1 1 -1 1 2 1 -1 2 8]
Candidate 3: Obj = 123.8613, schedule = [3 1 0 1 1 1 0 2 0 1 1 1 1 0 0 2 1 0 0 8]
Candidate 4: Obj = 127.0363, schedule = [2 1 0 1 0 2 1 1 0 1 0 1 2 0 0 2 1 0 1 8]
Best objective value found so far: 118.8710
Total points evaluated: 65
Iteration 3 completed in 3.93 seconds.
--- Iteration 4/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 2 0 1 0 1 2 -1 2 1 0 2 -1 2 0 2 0 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 1 0 1 0 2 1 0 1 0 1 1 0 1 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 1 0 1 0 1 1 1 0 2 1 1 0 0 2 0 2 -1 2 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 1 0 2 -1 1 2 0 0 2 1 0 2 -1 2 1 0 1 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 2 1 1 0 1 1 0 1 0 2 0 2 -1 1 9]
Best objective value found so far: 118.8710
Total points evaluated: 70
Iteration 4 completed in 4.20 seconds.
--- Iteration 5/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 1 1 1 1 1 1 1 0 0 2 1 0 1 0 9]
Candidate 1: Obj = 128.2031, schedule = [2 1 1 0 0 1 2 1 0 1 0 2 1 0 0 2 0 1 1 8]
Candidate 2: Obj = 125.0659, schedule = [1 1 1 1 0 1 2 0 1 1 1 1 0 1 0 2 0 1 1 8]
Candidate 3: Obj = 131.4179, schedule = [2 1 0 1 0 1 2 1 0 1 1 0 2 0 1 1 1 0 0 9]
Candidate 4: Obj = 132.0015, schedule = [2 0 1 1 0 1 1 1 1 1 1 1 0 0 2 1 1 0 0 9]
Best objective value found so far: 118.8710
Total points evaluated: 75
Iteration 5 completed in 4.52 seconds.
--- Iteration 6/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 0 1 1 0 2 0 1 0 1 2 1 1 -1 1 1 2 0 0 9]
Candidate 1: Obj = 123.5175, schedule = [2 1 1 1 0 2 0 1 0 1 1 2 0 0 2 1 0 0 1 8]
Candidate 2: Obj = 123.2617, schedule = [3 0 1 1 0 1 1 1 0 2 1 0 1 0 2 0 1 0 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 1 2 0 1 0 2 0 1 1 1 0 1 0 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 1 0 1 2 0 1 1 0 2 0 0 2 0 2 -1 1 9]
Best objective value found so far: 118.8710
Total points evaluated: 80
Iteration 6 completed in 4.32 seconds.
--- Iteration 7/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 2 0 0 2 0 2 -1 2 1 0 1 0 2 0 1 1 0 9]
Candidate 1: Obj = 133.3148, schedule = [2 1 1 0 0 2 0 1 0 1 2 0 1 0 2 0 1 0 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 2 1 -1 2 1 0 1 1 0 2 0 1 0 1 1 1 1 7]
Candidate 3: Obj = 124.1344, schedule = [3 1 1 0 1 1 1 0 0 2 0 1 1 1 0 2 0 0 2 7]
Candidate 4: Obj = 132.5415, schedule = [1 1 1 2 0 1 0 2 0 1 0 1 1 1 1 1 0 0 1 9]
Best objective value found so far: 118.8710
Total points evaluated: 85
Iteration 7 completed in 4.17 seconds.
--- Iteration 8/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 2 1 0 0 1 2 0 0 1 2 1 1 -1 2 1 0 0 2 7]
Candidate 1: Obj = 124.7751, schedule = [3 0 1 2 0 0 2 0 1 0 1 2 1 0 1 0 1 0 1 8]
Candidate 2: Obj = 122.3063, schedule = [2 1 2 0 0 1 2 1 0 1 1 1 1 0 1 0 1 0 2 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 2 1 1 0 0 1 1 0 1 1 2 0 0 1 1 2 -1 1 8]
Candidate 4: Obj = 132.8389, schedule = [2 0 2 1 0 1 0 1 1 0 1 2 0 0 1 2 0 0 2 8]
Best objective value found so far: 118.8710
Total points evaluated: 90
Iteration 8 completed in 4.06 seconds.
--- Iteration 9/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 2 0 1 1 0 2 -1 1 1 2 1 0 0 2 1 0 1 8]
Candidate 1: Obj = 132.6826, schedule = [2 1 0 2 0 0 2 1 0 0 1 1 1 1 0 2 0 0 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 1 2 1 0 1 0 1 2 0 0 1 1 1 1 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 1 1 0 2 0 2 0 0 1 1 2 -1 1 1 1 1 0 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 1 1 0 1 2 0 1 0 1 2 -1 2 0 1 0 2 7]
Best objective value found so far: 118.8710
Total points evaluated: 95
Iteration 9 completed in 3.84 seconds.
--- Iteration 10/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 0 2 1 0 0 2 1 -1 2 0 1 2 0 1 0 1 1 1 7]
Candidate 1: Obj = 131.3891, schedule = [2 0 1 1 1 1 1 0 0 2 0 1 2 0 1 0 1 1 0 9]
Candidate 2: Obj = 132.6630, schedule = [2 0 1 2 0 1 1 0 0 1 1 1 2 0 0 2 0 1 0 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 2 1 0 1 1 0 0 2 0 1 1 1 0 2 1 -1 2 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 2 1 1 0 1 1 0 1 1 0 1 1 0 1 9]
Best objective value found so far: 118.8710
Total points evaluated: 100
Iteration 10 completed in 4.38 seconds.
--- Iteration 11/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 119.2401, schedule = [2 1 2 0 1 1 1 0 0 1 2 1 0 0 1 1 1 1 1 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 1 2 0 0 2 1 0 0 1 2 1 1 -1 1 1 1 0 2 8]
Candidate 2: Obj = 129.2138, schedule = [1 1 2 0 0 2 1 0 0 1 2 1 0 0 1 1 1 1 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 2 0 0 2 0 2 0 0 1 2 1 -1 1 1 1 0 2 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 2 0 2 -1 1 1 2 0 0 2 0 1 0 2 7]
Best objective value found so far: 118.8710
Total points evaluated: 105
Iteration 11 completed in 4.10 seconds.
--- Iteration 12/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 1 2 0 0 2 1 0 1 0 1 1 1 0 2 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 1 2 0 2 -1 2 1 0 0 2 1 0 2 -1 2 0 1 1 0 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 1 0 1 0 2 1 -1 2 1 0 2 -1 2 0 1 1 0 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 1 1 1 -1 2 1 0 2 -1 1 2 0 0 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 0 2 -1 2 1 1 -1 2 1 0 1 0 2 1 0 0 2 8]
Best objective value found so far: 118.8710
Total points evaluated: 110
Iteration 12 completed in 4.41 seconds.
--- Iteration 13/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 118.2240, schedule = [3 0 1 2 0 0 1 1 1 0 2 0 1 0 2 0 1 1 1 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 1 1 2 -1 2 0 1 1 0 2 0 1 1 1 8]
Candidate 2: Obj = 117.9711, schedule = [3 0 1 2 0 0 1 1 1 0 1 1 1 0 2 0 2 0 1 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 2 0 1 2 -1 1 1 2 -1 2 0 1 1 1 1 0 1 0 2 8]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 2 0 0 1 1 1 1 0 2 0 0 2 0 2 -1 1 8]
Best objective value found so far: 117.9711
Total points evaluated: 115
Iteration 13 completed in 4.59 seconds.
--- Iteration 14/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 123.3827, schedule = [3 1 1 0 1 1 1 1 0 1 1 1 1 0 1 0 2 0 0 8]
Candidate 1: Obj = 132.5767, schedule = [2 1 1 0 1 0 2 1 0 0 2 1 1 0 1 1 1 0 0 9]
Candidate 2: Obj = 124.9912, schedule = [3 1 1 1 0 0 2 1 0 1 0 2 1 0 1 1 1 0 0 8]
Candidate 3: Obj = 119.2593, schedule = [3 1 1 1 0 1 1 0 1 1 1 0 2 0 0 2 1 0 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 0 1 2 0 0 2 1 1 -1 2 1 1 0 1 8]
Best objective value found so far: 117.9711
Total points evaluated: 120
Iteration 14 completed in 4.59 seconds.
--- Iteration 15/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 132.2421, schedule = [2 1 1 0 0 1 1 2 0 0 1 1 1 1 0 1 2 0 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 1 2 0 1 0 2 1 1 -1 1 2 0 1 1 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 0 2 0 1 0 1 1 0 1 1 2 -1 2 7]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 0 1 1 0 2 1 1 0 1 1 1 1 -1 2 0 2 -1 2 7]
Candidate 4: Obj = 117.2325, schedule = [3 0 1 1 0 2 0 1 1 1 1 0 2 0 0 2 1 0 1 7]
Best objective value found so far: 117.2325
Total points evaluated: 125
Iteration 15 completed in 4.36 seconds.
--- Iteration 16/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 1 0 0 2 1 0 1 1 0 1 2 0 0 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 0 1 0 2 1 0 2 0 0 1 2 0 0 8]
Candidate 2: Obj = 120.2212, schedule = [3 1 1 0 0 2 0 1 0 2 1 0 1 0 1 1 2 0 1 7]
Candidate 3: Obj = 122.3628, schedule = [2 2 0 1 0 2 0 1 0 2 1 0 1 1 0 1 2 0 0 8]
Candidate 4: Obj = 118.6656, schedule = [2 2 1 0 0 2 0 1 0 2 1 0 1 1 0 1 2 0 1 7]
Best objective value found so far: 117.2325
Total points evaluated: 130
Iteration 16 completed in 4.65 seconds.
--- Iteration 17/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 133.4606, schedule = [1 2 0 2 0 0 1 1 1 0 1 2 0 1 1 1 0 0 1 9]
Candidate 1: Obj = 124.3333, schedule = [3 1 0 2 0 0 2 0 1 1 0 2 0 0 2 1 0 1 0 8]
Candidate 2: Obj = 124.4456, schedule = [3 1 0 1 0 1 2 0 1 0 1 2 0 1 0 2 0 0 1 8]
Candidate 3: Obj = 122.1176, schedule = [2 1 1 2 0 0 1 1 0 2 0 2 0 1 1 1 0 0 2 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 2 1 1 2 -1 1 1 1 0 2 0 2 0 1 0 1 1 0 1 8]
Best objective value found so far: 117.2325
Total points evaluated: 135
Iteration 17 completed in 5.23 seconds.
--- Iteration 18/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 1 1 -1 1 2 0 1 1 0 2 0 0 1 2 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 2 0 0 2 1 0 1 1 0 1 1 0 1 2 1 -1 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 1 2 1 1 -1 1 2 0 1 0 1 2 0 0 1 2 1 0 0 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 1 1 1 1 0 2 0 1 0 1 2 0 0 1 2 1 -1 2 8]
Candidate 4: Obj = 132.6467, schedule = [2 0 2 1 0 1 1 0 1 1 0 2 0 0 1 2 0 0 2 8]
Best objective value found so far: 117.2325
Total points evaluated: 140
Iteration 18 completed in 5.33 seconds.
--- Iteration 19/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 132.7635, schedule = [1 2 1 0 1 1 0 1 1 0 1 2 0 1 0 1 1 0 1 9]
Candidate 1: Obj = 134.0218, schedule = [1 2 0 2 0 1 0 1 1 0 2 0 1 0 1 2 0 0 1 9]
Candidate 2: Obj = 132.6921, schedule = [1 2 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 0 1 9]
Candidate 3: Obj = 134.3468, schedule = [1 2 1 1 0 1 0 1 1 1 0 2 0 0 1 2 0 0 1 9]
Candidate 4: Obj = 131.9700, schedule = [2 1 1 1 0 1 0 1 0 1 1 2 0 1 0 1 1 1 0 9]
Best objective value found so far: 117.2325
Total points evaluated: 145
Iteration 19 completed in 5.02 seconds.
--- Iteration 20/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 124.8142, schedule = [2 1 0 2 0 1 0 2 0 0 1 2 1 0 1 1 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 1 1 0 1 1 0 1 2 -1 2 0 1 1 0 9]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 2 0 1 1 0 0 2 0 1 2 0 0 2 1 -1 1 9]
Candidate 3: Obj = 133.0355, schedule = [2 1 0 2 0 1 1 1 0 0 1 2 0 0 1 2 0 1 0 9]
Candidate 4: Obj = 127.4777, schedule = [2 1 1 1 0 0 2 0 0 1 1 1 1 1 0 2 1 0 1 8]
Best objective value found so far: 117.2325
Total points evaluated: 150
Iteration 20 completed in 5.28 seconds.
--- Iteration 21/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 3 0 2 1 -1 1 2 1 0 0 1 1 1 1 0 1 2 -1 2 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 0 2 0 1 0 1 1 1 0 2 0 0 2 7]
Candidate 2: Obj = 118.8553, schedule = [2 2 1 0 0 2 0 1 1 0 1 2 0 0 2 1 0 1 1 7]
Candidate 3: Obj = 123.2719, schedule = [3 1 0 1 1 0 2 0 1 0 1 1 1 0 1 1 1 0 1 8]
Candidate 4: Obj = 122.1658, schedule = [3 1 0 1 1 0 2 1 0 1 1 0 1 1 0 1 1 1 0 8]
Best objective value found so far: 117.2325
Total points evaluated: 155
Iteration 21 completed in 5.98 seconds.
--- Iteration 22/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 121.0466, schedule = [3 0 1 1 1 1 0 1 1 0 2 0 1 1 1 0 1 0 2 7]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 1 1 -1 1 2 1 0 0 2 1 1 -1 1 1 1 0 2 8]
Candidate 2: Obj = 10000000000.0000, schedule = [ 3 0 1 2 0 1 0 2 -1 1 1 2 1 0 1 1 1 0 1 7]
Candidate 3: Obj = 134.6840, schedule = [1 2 1 0 0 1 1 2 0 1 0 1 2 0 0 2 0 0 1 9]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 2 -1 2 1 1 0 1 0 2 0 1 0 1 1 0 1 8]
Best objective value found so far: 117.2325
Total points evaluated: 160
Iteration 22 completed in 5.36 seconds.
--- Iteration 23/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 2 1 1 1 0 0 1 2 -1 1 2 0 1 1 0 2 1 0 1 8]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 0 2 0 1 1 0 1 0 2 0 2 0 1 1 0 2 -1 2 8]
Candidate 2: Obj = 128.1131, schedule = [2 1 0 1 0 1 2 1 0 0 2 0 2 0 0 2 0 1 1 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 1 1 0 1 0 1 2 -1 1 1 1 2 -1 1 1 1 1 1 7]
Candidate 4: Obj = 10000000000.0000, schedule = [ 3 0 1 1 1 1 1 1 -1 1 2 1 0 1 0 2 1 0 1 7]
Best objective value found so far: 117.2325
Total points evaluated: 165
Iteration 23 completed in 5.23 seconds.
--- Iteration 24/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 1 2 0 0 1 2 0 1 0 2 1 1 0 0 9]
Candidate 1: Obj = 10000000000.0000, schedule = [ 2 1 1 2 -1 1 2 1 -1 2 1 0 1 0 1 1 2 0 1 7]
Candidate 2: Obj = 10000000000.0000, schedule = [ 2 1 0 1 0 1 2 1 0 1 1 0 2 0 1 0 2 -1 1 9]
Candidate 3: Obj = 10000000000.0000, schedule = [ 1 2 0 2 0 0 2 1 -1 2 1 1 1 0 0 1 2 -1 1 9]
Candidate 4: Obj = 10000000000.0000, schedule = [ 1 1 1 2 -1 1 1 1 0 2 1 0 1 0 2 0 2 0 1 8]
Best objective value found so far: 117.2325
Total points evaluated: 170
Iteration 24 completed in 5.52 seconds.
--- Iteration 25/25 ---
Fitting GP model...
GP model fitted.
Optimizing acquisition function (LCB)...
Selected 5 candidate(s).
Candidate 0: Obj = 10000000000.0000, schedule = [ 1 2 1 1 0 1 0 2 0 0 2 0 1 0 2 1 1 -1 2 8]
Candidate 1: Obj = 124.5451, schedule = [3 0 2 0 0 1 2 0 1 0 2 0 2 0 1 0 1 0 1 8]
Candidate 2: Obj = 124.2147, schedule = [3 0 1 1 1 0 2 0 0 2 0 2 1 0 0 2 1 0 0 8]
Candidate 3: Obj = 10000000000.0000, schedule = [ 3 1 1 1 -1 2 0 1 0 2 0 2 0 0 1 1 1 1 1 7]
Candidate 4: Obj = 123.0855, schedule = [2 2 0 1 1 1 1 0 0 1 2 0 1 0 1 1 2 0 0 8]
Best objective value found so far: 117.2325
Total points evaluated: 175
Iteration 25 completed in 6.93 seconds.
--- Optimization Finished ---
Total evaluations: 175
Best Objective Value Found: 117.23250201869716
Best U vector Found: [0 1 0 0 0 0 1 0 0 1 1 1 0 1 1 0 1 1 1 1]
--- Verification ---
Is the best U feasible? True
Resulting Y vector for best U: [3 0 1 1 0 2 0 1 1 1 1 0 2 0 0 2 1 0 1 7]
Objective value (recalculated): 117.23256.5 Results
The initial schedule, derived using the Bailey-Welch method (bailey1952study), serves as a baseline. The objective function \(C(\mathbf{x})\) combines Expected Waiting Time (\(EWT\)) and Expected Staff Penalty (\(ESP\)). Lower values of \(C(\mathbf{x})\) are preferable. Each experiment consisted of \(N_{INITIAL} = 20\) initial random evaluations followed by \(N_{ITERATIONS} =20\) Bayesian optimization iterations, with \(BATCH\_SIZE_q = 5\) evaluations per iteration, totaling approximately \(20 + 20 \times 5 = 120\) evaluations per experiment. The optimization operates on the binary perturbation vector \(\mathbf{U}\), using the HED embedding (Deshwal et al. 2023).
The key performance metric is the best (minimum) objective function value found.
6.5.1 Experiment 1: CBO with Expected Improvement (EI)
- Initial Best Objective (after random search): 118.87100166955729
- Final Best Objective Found: 114.96364292044709
- Total Evaluations: 175
- Best Perturbation Vector \(\mathbf{U}_{EI}^*\): array([0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1])
- Resulting Optimal Schedule \(\mathbf{x}_{EI}^*\): array([2, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 2, 0, 1, 7])
6.5.2 Experiment 2: CBO with Lower Confidence Bound (LCB) - \(\kappa =\) 2.576
- Initial Best Objective (after random search): 118.87100166955729
- Final Best Objective Found: 117.23250201869716
- Total Evaluations: 175
- Best Perturbation Vector \(\mathbf{U}_{LCB\_fixed}^*\): array([0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1])
- Resulting Optimal Schedule \(\mathbf{x}_{LCB\_fixed}^*\): array([3, 0, 1, 1, 0, 2, 0, 1, 1, 1, 1, 0, 2, 0, 0, 2, 1, 0, 1, 7])
6.5.3 Summary of Best Objectives
| Experiment | Best Objective \(C(\mathbf{x}^*)\) |
|---|---|
| CBO with EI | 114.96364292044709 |
| CBO with LCB (Fixed \(\kappa=\) 2.576) | 117.23250201869716 |
6.6 Discussion
The experiments aimed to compare two CBO strategies, leveraging the HED embedding technique (Deshwal et al. 2023), for optimizing the outpatient appointment scheduling problem formulated by Kaandorp and Koole (2007). Both methods successfully improved upon their respective initial random search results, demonstrating the applicability of BO with HED to this combinatorial problem.
Performance Comparison:
Hypothesis Evaluation:
- Hypothesis 1
- Hypothesis 2
Exploration vs. Exploitation:
Computational Effort:
Limitations and Future Work:
- The optimality guarantee mentioned by Kaandorp and Koole (2007) applies to their specific local search algorithm operating directly on the schedule space \(\mathcal{F}\), leveraging multimodularity. Our BO approach operates on the perturbation vector space \(\mathbf{U}\) via HED embeddings. While BO aims for global optimization, it doesn’t inherit the same theoretical guarantee of finding the global optimum as the original local search, especially given the stochastic nature of GP modeling and acquisition function optimization.
- The performance is likely sensitive to BO hyperparameters (dictionary size \(m\), \(\kappa\) values, number of candidates for acquisition optimization).
- Further investigation into different dictionary construction methods (e.g., binary wavelets as mentioned in Deshwal et al., 2023) or adaptive \(\kappa\) schedules could be beneficial.
In conclusion, applying CBO with HED embeddings appears promising for this scheduling problem. The LCB acquisition function with a fixed, well-chosen \(\kappa\) demonstrated the best performance in this study.
6.7 Timeline
- Experiment Setup and Code Implementation: 30-04-2025
- Results Analysis and Report Compilation: 07-05-2025
6.8 References
Aglietti, Virginia, Ira Ktena, Jessica Schrouff, Eleni Sgouritsa,
Francisco J. R. Ruiz, Alan Malek, Alexis Bellot, and Silvia Chiappa.
2024. “FunBO: Discovering
Acquisition Functions for
Bayesian Optimization with
FunSearch.” arXiv. https://doi.org/10.48550/arXiv.2406.04824.
Altman, Eitan, Bruno Gaujal, and Arie Hordijk. 2000.
“Multimodularity, Convexity, and
Optimization Properties.”
Mathematics of Operations Research 25 (2): 324–47. https://www.jstor.org/stable/3690584.
Bailey, Norman TJ. 1952. “A Study of Queues and Appointment
Systems in Hospital Out-Patient Departments, with Special Reference to
Waiting-Times.” Journal of the Royal Statistical Society
Series B: Statistical Methodology 14 (2): 185–99.
Balandat, Maximilian, Brian Karrer, Daniel R. Jiang, Samuel Daulton,
Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. 2020.
“BoTorch: A Framework for Efficient
Monte-Carlo Bayesian Optimization.” In Advances in
Neural Information Processing Systems 33. http://arxiv.org/abs/1910.06403.
Deshwal, Aryan, Sebastian Ament, Maximilian Balandat, Eytan Bakshy,
Janardhan Rao Doppa, and David Eriksson. 2023. “Bayesian
Optimization over
High-Dimensional Combinatorial
Spaces via Dictionary-Based
Embeddings.” In Proceedings of The
26th International Conference on
Artificial Intelligence and
Statistics, 7021–39. PMLR. https://proceedings.mlr.press/v206/deshwal23a.html.
Gonzalez, Javier, Zhenwen Dai, Andreas Damianou, and Neil D. Lawrence.
2017. “Preferential Bayesian
Optimization.” arXiv. https://doi.org/10.48550/arXiv.1704.03651.
Kaandorp, Guido C., and Ger Koole. 2007. “Optimal Outpatient
Appointment Scheduling.” Health Care Management Science
10 (3): 217–29. https://doi.org/10.1007/s10729-007-9015-x.
Swersky, Kevin, Yulia Rubanova, David Dohan, and Kevin Murphy. 2020.
“Amortized Bayesian Optimization over
Discrete Spaces.” In Proceedings of
the 36th Conference on Uncertainty in
Artificial Intelligence
(UAI), 769–78. PMLR. https://proceedings.mlr.press/v124/swersky20a.html.
Zhang, Xinyu, Daolang Huang, Samuel Kaski, and Julien Martinelli. 2025.
“PABBO: Preferential
Amortized Black-Box
Optimization.” arXiv. https://doi.org/10.48550/arXiv.2503.00924.