Impact Statement

Providing a practical method for the optimization of expensive experiments and simulations under changing and externally given environmental conditions allows solving engineering problems that better mimic the challenges found in real-world applications. It makes results robust when transferring them from the lab to the real world as it removes simplistic assumptions that fail to hold in realistic scenarios. Instead of running one experiment multiple times under different environmental conditions, ENVBO, the proposed method, requires only a single optimization run that shares all available information about the controllable variables and the environmental conditions. This decreases the number of observations needed and lowers the cost of optimization, which is the main challenge of expensive black-box optimization.

2. Related work

The methodology in this article is in the area of Bayesian optimization, whose origin can be traced back to the early 1960s when Kushner (Reference Kushner1962) and Kushner (Reference Kushner1964) introduced an optimization algorithm on a one-dimensional noisy function using a Gaussian process as a surrogate model and sampling heuristics that later evolved into the upper confidence bound (UCB) and probability of improvement acquisition functions. Research from 1971 onward linked Bayesian optimization predominantly to the optimization of expensive multimodal functions and introduced new heuristics—notably expected improvement (EI; Šaltenis Reference Šaltenis1971) and knowledge gradient (Močkus Reference Močkus1972, Reference Močkus1974, Reference Močkus1989). Limitations in computational capabilities at the time resulted in using cheaper models, such as Wiener and Ornstein–Uhlenbeck processes, for which both heuristics collapse onto each other (Garnett Reference Garnett2023). The evolution of technology allowed computationally expensive models such as Gaussian processes to be used in the late 1990s, which incentivized reevaluation of the methods proposed two decades earlier. Schonlau (Reference Schonlau1997) and Jones et al. (Reference Jones, Schonlau and Welch1998) introduced their efficient global optimization (EGO) algorithm and studied EI in a deterministic setting. Auer (Reference Auer2003) and Srinivas et al. (Reference Srinivas, Krause, Kakade and Seegre2010) reexamined UCB and studied its theoretical regret for the first time, while Frazier and Powell (Reference Frazier and Powell2007) and Scott et al. (Reference Scott, Frazier and Powell2011) studied the knowledge gradient for discrete and continuous input spaces, respectively. Villemonteix et al. (Reference Villemonteix, Vázquez and Walter2006) introduced an information-based strategy into the Bayesian optimization framework that was later studied and further developed into entropy search (Hennig and Schuler Reference Hennig and Schuler2011), predictive entropy search (Hernández-Lobato et al. Reference Hernández-Lobato, Hoffman and Ghahramani2014), and max-value entropy search (Wang and Jegelka Reference Wang and Jegelka2017). Comprehensive surveys are available from Brochu et al. (Reference Brochu, Cora and de Freitas2010), Shahriari et al. (Reference Shahriari, Swersky, Wang, Adams and De Freitas2015) and Frazier (Reference Frazier2018). Closely related to Bayesian optimization is the work on multiarmed bandits investigating a similar problem but in a discrete parameter space (Garnett Reference Garnett2023; Berry and Fristedt Reference Berry and Fristedt1985; Lattimore and Szepesvári Reference Lattimore and Szepesvári2020).

Bayesian optimization has been applied to many scientific fields, predominantly in the science, technology, engineering, and mathematics areas. Applications include fields such as material science (Frazier and Wang Reference Frazier and Wang2015; Packwood Reference Packwood2017) and drug discovery (Negoescu et al. Reference Negoescu, Frazier and Powell2011) in chemistry, physics (Duris et al. Reference Duris, Kennedy, Hanuka, Shtalenkova, Edelen, Egger, Cope and Ratner2019), and biology (Li et al. Reference Li, Reyes, Vazquez-Anderson, Wang, Contreras and Powell2018). In engineering, it has been used in civil engineering (Gramacy et al. Reference Gramacy, Gray, Digabel, Lee, Ranjan, Wells and Wild2016; Garnett et al. Reference Garnett, Osborne and Roberts2010), electrical engineering (Lyu et al. Reference Lyu, Xue, Yang, Yan, Hong, Zeng and Zhou2018), mechanical engineering (Sterling et al. Reference Sterling, Sterling, Zhang and Chen2015), and robotics (Martinez-Cantin et al. Reference Martinez-Cantin, de Freitas, Brochu, Castellanos and Doucet2009; Calandra et al. Reference Calandra, Seyfarth, Peters and Deisenroth2015). Furthermore, Bayesian optimization is often used in machine learning and artificial intelligence, for example, to configure algorithms (Hutter et al. Reference Hutter, Hoos and Leyton-Brown2011) and to tune hyperparameters (Snoek et al. Reference Snoek, Larochelle and Adams2012). Garnett (Reference Garnett2023) provides a thorough list of applications.

This article focuses on a particular case of Bayesian optimization (Frazier Reference Frazier2018) that is referred to as random environmental conditions (Chang et al. Reference Chang, Williams, Bhalla, Belknap, Santner, Notz and Bartel2001), multitask optimization (Swersky et al. Reference Swersky, Snoek and Adams2013), and contextual optimization (Krause and Ong Reference Krause and Ong2011). Methods differ by the type of input parameters (discrete or continuous), the number of allowed contexts (finite or infinite), and the overall goal of the optimization (one optimal solution for all contexts or one optimal solution for each context). When the number of contexts is infinite, finding one optimal solution for each context corresponds to finding a function that returns optimal inputs based on the context. This is the main objective of this article. Gaussian process-based methods aiming to optimize problems with different contexts—that is, with controllable and environmental variables—considered in the past can be classified primarily into two types based on their optimization goal.

The first type aims to find one solution that yields the best result for all tasks/contexts. It is assumed that environmental variables take values according to a probability distribution. For example, if the environment is defined by its temperature, the current temperature will follow a distribution. Optimization aims to find the optimizer that maximizes the objective function, taking into account the likelihood of the different environments. Different approaches assume the distribution of the environmental variables to be discrete (Williams et al. Reference Williams, Santner and Notz2000; Swersky et al. Reference Swersky, Snoek and Adams2013), continuous (Groot et al. Reference Groot and Birlutiu2010; Xie et al. Reference Xie, Frazier, Sankaran, Marsden and Elmohamed2012), or both (Toscano-Palmerin and Frazier Reference Toscano-Palmerin and Frazier2018). These methods assume that parameters and contexts can be selected during optimization, and many use a sequential approach, where parameters for the next candidate are selected before the context is chosen. Only Toscano-Palmerin and Frazier (Reference Toscano-Palmerin and Frazier2018) present a method that chooses parameters and contexts jointly and can optimize problems where contexts are uncontrollable and given randomly. Chang et al. (Reference Chang, Williams, Santner, Notz and Bartel1999) and Chang et al. (Reference Chang, Williams, Bhalla, Belknap, Santner, Notz and Bartel2001) use this type of method to design a femoral component for hip replacements conditional on joint force orientation and cancellous bone properties. The aim is to find one optimal design for a wide demographic with varying characteristics.

The second type of method aims to find one solution for each context. Thus, interest lies not in finding one global optimum but multiple optima, one for each combination of environmental conditions. Pearce and Branke (Reference Pearce and Branke2017), Char et al. (Reference Char, Chung, Neiswanger, Kandasamy, Nelson, Boyer, Kolemen and Schneider2019), and Chung et al. (Reference Chung, Char, Neiswanger, Kandasamy, Nelson, Boyer, Kolemen and Schneider2020) consider multiple discrete tasks, while Ginsbourger et al. (Reference Ginsbourger, Baccou, Chevalier, Perales, Garland and Monerie2014) and Pearce and Branke (Reference Pearce and Branke2018) consider continuous environmental variables. These methods are closely related to the objective of this article. However, they have one important distinction. While the environmental values are given externally in real-world applications, it is assumed that they can be set to any desired values in the experiments and simulations above. This deviates from our problem formulation, where we explicitly regard problems with uncontrollable environmental variables—in experiments, simulations, and the real world.

The research in this article is closely related to Krause and Ong (Reference Krause and Ong2011), who modified UCB (Srinivas et al. Reference Srinivas, Krause, Kakade and Seegre2010) to be suited to optimization with externally given environmental conditions and derived theoretical bounds for its contextual regret. In contrast to Krause and Ong (Reference Krause and Ong2011), this article considers improvement-based acquisition functions, that is, EI (Jones et al. Reference Jones, Schonlau and Welch1998) and log EI (LogEI; Ament et al. Reference Ament, Daulton, Eriksson, Balandat and Bakshy2023), gives a detailed description of the practical implementation of Bayesian optimization with environmental conditions and provides all code at https://github.com/mikediessner/environmental-conditions-BO. In addition, the reported approach makes fewer assumptions than Krause and Ong (Reference Krause and Ong2011), who focus on a linear and additive covariance structure for the environmental variables. These assumptions can be restrictive as no information about the objective function might be available—as is often the case with expensive, derivative-free black-box optimization. Thus, a widely applicable strategy is required that this article aims to provide.

The area of multitask optimization is also related to multi-objective optimization. In multi-objective optimization, the aim is no longer to optimize one single objective function but rather multiple objective functions simultaneously. A global optimum for all objective functions is only possible if it happens to be at the same values for all inputs for all functions. Thus, the main goal is to find the best trade-off between the objective functions (Garnett Reference Garnett2023; Frazier Reference Frazier2018). This is related to the first problem above, which aims to find a trade-off between all contexts. However, objective functions in multi-objective optimization are generally assumed not to be correlated, while correlations are leveraged for contexts. Acquisition functions such as EI (Emmerich et al. Reference Emmerich, Giannakoglou and Naujoks2006; Ponweiser et al. Reference Ponweiser, Wagner, Biermann and Vincze2008; Yang et al. Reference Yang, Emmerich, Deutz and Bäck2019b; Yang et al. Reference Yang, Emmerich, Deutz and Bäck2019a), predictive entropy search (Hernández-Lobato et al. Reference Hernández-Lobato, Gelbart, Hoffman, Adams and Ghahramani2015) and max-value entropy search (Belakaria et al. Reference Belakaria, Deshwal and Doppa2020; Fernández-Sánchez et al. Reference Fernández-Sánchez, Garrido-Merchán and Hernández-Lobato2023) have been extended to allow multi-objective optimization.

3. Methodology

This section provides a brief overview of the fundamentals of Bayesian optimization, particularly the surrogate modeling via Gaussian processes and the acquisition functions used to guide the optimization. It then introduces the novel conditional Bayesian optimization algorithm that enables optimization with uncontrollable environmental variables.

3.1. Bayesian optimization

Consider the d-dimensional maximization problem

(1) $$ {\boldsymbol{x}}^{\ast }=\underset{\boldsymbol{x}\in \mathcal{X}}{\mathrm{argmax}}\hskip0.3em f\left(\boldsymbol{x}\right), $$

where $ \mathcal{X} $ is a continuous input space that is bounded by a hyper-rectangle such that $ \mathcal{X}\in {\left[a,b\right]}^d $ with $ a,b\in \mathrm{\mathbb{R}} $ . The objective function $ f\left(\mathbf{x}\right) $ typically has three properties. First, it is a black-box function that can be provided with an input vector $ {\mathbf{x}}_i $ and allows the observation of the scalar output $ {y}_i $ . Beyond this, no other information can be inferred from the function. Second, the function is expensive to evaluate, and significant time, resources, and money costs are generated at each evaluation. Third, the function typically lacks a derivative or is too expensive to compute (Frazier Reference Frazier2018). Any noise $ \unicode{x025B} $ introduced during the evaluation of the objective function $ f\left(\mathbf{x}\right) $ is assumed to be independent and identically distributed Gaussian noise $ \unicode{x025B} \sim \mathcal{N}\left(0,{\sigma}^2\right) $ such that an observation can be defined as $ {y}_i=f\left({\mathbf{x}}_i\right)+\unicode{x025B} $ . Multiple observation pairs consisting of inputs and outputs are defined as $ {\mathcal{D}}_n={\left\{\left({\mathbf{x}}_i,{y}_i\right)\right\}}_{i=1}^n $ . In this article, $ {\mathbf{X}}_n={\left\{{\mathbf{x}}_i\right\}}_{i=1}^n $ and $ {\mathbf{y}}_n={\left\{{y}_i\right\}}_{i=1}^n $ are used to describe all training inputs and their corresponding outputs.

Bayesian optimization (Močkus Reference Močkus1974, Reference Močkus1989; Žilinskas Reference Žilinskas1975Jones et al. Reference Jones, Schonlau and Welch1998; Srinivas et al. Reference Srinivas, Krause, Kakade and Seegre2010; Snoek et al. Reference Snoek, Larochelle and Adams2012; Shahriari et al. Reference Shahriari, Swersky, Wang, Adams and De Freitas2015; Frazier Reference Frazier2018; Gramacy Reference Gramacy2020) is an optimization algorithm based on surrogate modeling to solve expensive problem (1) in a minimum number of function evaluations. The expensive and opaque nature of these problems requires a sample-efficient optimization algorithm that keeps costs low to make the optimization feasible. Bayesian optimization has emerged as a prime candidate using a surrogate model $ \mathcal{M} $ to represent the unknown objective function $ f\left(\mathbf{x}\right) $ . The mean of the surrogate model’s predictive distribution is then used to compute an acquisition criterion $ \alpha \left(\cdot \right) $ that guides the optimization process by proposing new input points to be evaluated by the objective function. Bayesian optimization is a sequential optimization algorithm performed in a feedback loop, as illustrated in Algorithm 1. This loop consists of three steps. First, the surrogate model is fitted to the training data $ {\mathcal{D}}_n $ . Second, the next candidate point $ {\mathbf{x}}_{n+1} $ is computed by maximizing the acquisition criterion $ \alpha $ . Third, the objective function evaluates the new candidate point, and its output $ {y}_{n+1} $ is observed. The loop then starts again, adding this newly observed candidate point to the training set. Thus, Bayesian optimization gathers more information sequentially with each loop. The algorithm stops when a predefined evaluation budget $ N $ is exhausted and returns the data pair with the highest observation as its solution.

Algorithm 1 Basic Bayesian optimization algorithm.

Require: Evaluation budget $ N $ , number of initial points $ {n}_0 $ , surrogate model $ \mathcal{M} $ , and acquisition function $ \alpha $ .

Sample $ {n}_0 $ initial training data points $ {\mathbf{X}}_0 $ via a space-filling design (McKay et al. Reference McKay, Beckman and Conover1979) and gather observations $ {\mathbf{y}}_0 $ .

Set $ n=0 $ .

while $ n\le N-{n}_0 $ do.

Fit surrogate model $ \mathcal{M} $ to training data $ {\mathcal{D}}_n=\left\{{\mathbf{X}}_n,{\mathbf{y}}_n\right\} $ .

Find $ {\mathbf{x}}_{n+1} $ that maximizes an acquisition criterion $ \alpha $ based on model $ \mathcal{M} $ , that is, solve $ {\mathrm{argmax}}_{\mathbf{x}}\alpha \left(\mathbf{x}\right) $ .

Evaluate $ {\mathbf{x}}_{n+1} $ , observing $ {y}_{n+1} $ .

Increment $ n $ .

end while.

return Point $ {\mathbf{x}}^{\ast } $ with the highest observation $ {y}^{\ast } $ .

Algorithm 1 is illustrated in Figure 1. The objective function (dashed line) is a simple one-dimensional function with one local optimum and one global optimum at $ x=8 $ . The algorithm is initialized with three observations (dark blue dots), and the surrogate model is fitted, providing its prediction (red line) and the corresponding uncertainty (blue area) in the form of 95% confidence intervals around the prediction. This model is used to compute the acquisition criterion (orange area)—in this case, EI (see Section 3.1.2)—that, when maximized, provides the next candidate point (dashed red line) to be evaluated from the objective function. Bayesian optimization is run for eight iterations, and the surrogate model is updated with each newly evaluated candidate point until the algorithm finds the optimum at iteration six.

Figure 1. Bayesian optimization is applied to a one-dimensional function with one local and one global maximum. EI is used as the acquisition function. The input space is bounded by $ \left[0,10\right] $ .

3.1.1. Gaussian process

Different methods, such as random forest regression (Hutter et al. Reference Hutter, Hoos and Leyton-Brown2011; Lindauer et al. Reference Lindauer, Eggensperger, Feurer, Biedenkapp, Deng, Benjamins, Ruhkopf, Sass and Hutter2022), density ratio estimation (Bergstra et al. Reference Bergstra, Bardenet, Bengio and Kégl2011; Tiao et al. Reference Tiao, Klein, Seeger, Bonilla, Archambeau and Ramos2021), and neural networks (Snoek et al. Reference Snoek, Rippel, Swersky, Kiros, Satish, Sundaram, Patwary and Prabhat Adams2015; Springenberg et al. Reference Springenberg, Klein, Falkner and Hutter2016; Maraval et al. Reference Maraval, Zimmer, Grosnit, Tutunov, Wang and Ammar2022) have been considered for surrogate modeling in Bayesian optimization. However, Gaussian process regression (Rasmussen and Williams Reference Rasmussen and Williams2006; Gramacy Reference Gramacy2020) has seen the widest adaption for the surrogate model $ \mathcal{M} $ and has been researched extensively since the 1970s (Brochu et al. Reference Brochu, Cora and de Freitas2010). The four main reasons for this adaption are: first, Gaussian processes are a flexible modeling strategy, capable of representing a wide range of possible continuous objective functions $ f\left(\mathbf{x}\right) $ (Brochu et al. Reference Brochu, Cora and de Freitas2010; Frazier Reference Frazier2018; Gramacy Reference Gramacy2020; Garnett Reference Garnett2023). Second, they allow incorporating prior knowledge about the objective function by choosing the prior mean function and the prior covariance kernel (Brochu et al. Reference Brochu, Cora and de Freitas2010; Frazier Reference Frazier2018; Gramacy Reference Gramacy2020). The latter is particularly important as it controls the smoothness of the possible representations of the objective function. Third, they return a prediction with corresponding uncertainty, both required for computing acquisition functions (Garnett Reference Garnett2023; Gramacy Reference Gramacy2020). Fourth, they are appropriate for continuous objective functions, which are the main focus of Bayesian optimization, and work well with the limited amounts of data available in expensive black-box optimization (see Section 3.1). Due to these advantages, Gaussian processes are typically used in Bayesian optimization; thus, this article considers Bayesian optimization based on Gaussian process regression. Rasmussen and Williams (Reference Rasmussen and Williams2006) provide a comprehensive overview of Gaussian process regression and its implementation.

A Gaussian process is a nonparametric regression model that returns a prediction and its corresponding uncertainty for an unobserved point and can be described as a distribution of possible objective functions. Mathematically, it is a finite collection of random variables that have a joint Gaussian distribution. A Gaussian process only requires a prior mean function $ {\mu}_0\left(\mathbf{x}\right):\mathcal{X}\hskip0.3em \mapsto \hskip0.3em \mathrm{\mathbb{R}} $ and a prior covariance kernel $ {\Sigma}_0\left(\mathbf{x},{\mathbf{x}}^{\prime}\right):\mathcal{X}\times \mathcal{X}\hskip0.35em \mapsto \hskip0.35em \mathrm{\mathbb{R}} $ to be fully defined. The mean vector $ m\left({\mathbf{X}}_n\right)\hskip0.3em := \hskip0.4em {\mu}_0\left({\mathbf{X}}_n\right) $ and the $ n\times n $ covariance matrix $ K\left({\mathbf{X}}_n,{\mathbf{X}}_n\right)\hskip0.3em := \hskip0.4em {\Sigma}_0\left({\mathbf{X}}_n,{\mathbf{X}}_n\right) $ specify the multivariate normal distribution, also called the prior distribution, as

(2) $$ f\left({\mathbf{X}}_n\right)\sim \mathcal{N}\left(m\left({\mathbf{X}}_n\right),K\left({\mathbf{X}}_n,{\mathbf{X}}_n\right)\right). $$

This study follows Snoek et al. (Reference Snoek, Larochelle and Adams2012) and chooses the constant mean function given in Equation (3) as the prior mean function $ {\mu}_0\left(\cdot \right) $ and the Matérn 5/2 kernel presented in Equation (4) as the prior covariance function $ {\Sigma}_0\left(\cdot, \cdot \right) $ . Snoek et al. (Reference Snoek, Larochelle and Adams2012) argue that the Matérn 5/2 kernel is better suited for the optimization of realistic black-box problems than the often-used radial basis function kernel (Equation (5)) because the latter assumes the objective function to be very smooth. This assumption might be too strong for many problems, making the Matérn kernel a good alternative. However, the main advantage of using Gaussian processes as the surrogate model is the ability to incorporate prior knowledge about the objective function. Thus, the choice of the covariance kernel should always be informed by the available information about the given problem. For cases where information about the smoothness of the objective function is unavailable—which might be common for black-box problems—, Snoek et al. (Reference Snoek, Larochelle and Adams2012) make a case in favor of the Matérn 5/2 kernel over the radial basis function kernel as a default.

The Matérn 5/2 kernel uses the distance between inputs $ r=\mid \mathbf{x}-{\mathbf{x}}^{\prime}\mid $ to compute the uncertainty around its prediction. The kernel uses the output scale $ {\sigma}_f^2 $ to scale the covariance, where smaller values correspond to a smaller deviation from its mean. The characteristic length scale $ l $ quantifies the extent to which outputs are correlated when moving along an input axis. Smaller length scales $ l $ mean shorter correlation lengths and more variable functions, while larger length scales $ l $ correspond to longer correlation lengths and more constant functions (Gramacy Reference Gramacy2020; Rasmussen and Williams Reference Rasmussen and Williams2006).

(3) $$ \mu \left(\mathbf{x}\right)=c $$

(4) $$ {\Sigma}_{\mathrm{Mat}{\mathrm{e}}^{\prime}\mathrm{rn}}\left(\mathbf{x},{\mathbf{x}}^{\prime}\right)={\sigma}_f^2\left(1+\frac{\sqrt{5}r}{l}+\frac{5{r}^2}{3{l}^2}\right)\exp \left(-\frac{\sqrt{5}r}{l}\right) $$

The radial basis function kernel given in Equation (5) is another popular alternative for the covariance function. It is much smoother than the Matérn kernel and thus is only suited for Bayesian optimization when it can be assumed that the underlying objective function is necessarily smooth (Snoek et al. Reference Snoek, Larochelle and Adams2012).

(5) $$ {\Sigma}_{\mathrm{RBF}}\left(\mathbf{x},{\mathbf{x}}^{\prime}\right)={\sigma}_f^2\exp \left(-\frac{r^2}{2{l}^2}\right) $$

Covariance kernels can use one characteristic length scale $ l $ for all input dimensions $ d $ or $ d $ length scales $ \mathbf{l}={\left\{{l}_i\right\}}_{i=1}^d $ , one for each input dimension $ d $ . While the former has less computational overhead due to only requiring one length scale for all dimensions, it is less flexible as it assigns the same correlation to each dimension. The latter considers each dimension individually and assigns an individual length scale to each dimension. This is known as automatic relevance determination (Neal Reference Neal1996) as the Gaussian process assigns larger length scales to more constant dimensions, making them less influential in the computation of the covariance matrix. In contrast, inputs that are variable and change quickly are assigned smaller length scales, increasing their importance when computing the covariance kernel, as changes in these inputs generally affect the prediction significantly. Generally, the inverse of the length scales $ \mathbf{l} $ indicates the relevance of the corresponding input (Rasmussen and Williams Reference Rasmussen and Williams2006).

The Gaussian process described above has some hyperparameters $ \theta $ that can be estimated from the training data by maximizing the log-marginal likelihood in Equation (6) via maximum likelihood estimation (Rasmussen and Williams Reference Rasmussen and Williams2006).

(6) $$ {\displaystyle \begin{array}{l}\log p\left({\mathbf{y}}_n|{\mathbf{X}}_n\right)=-\frac{1}{2}{\left({\mathbf{y}}_n-m\left({\mathbf{X}}_n\right)\right)}^{\top }{\left[K\left({\mathbf{X}}_n,{\mathbf{X}}_n\right)+{\sigma}_y^2I\right]}^{-1}\left({\mathbf{y}}_n-m\left({\mathbf{X}}_n\right)\right)\\ {}\hskip10.1em -\frac{1}{2}\log \mid K\left({\mathbf{X}}_n,{\mathbf{X}}_n\right)+{\sigma}_y^2I\mid -\frac{n}{2}\log 2\pi \end{array}} $$

Besides the constant $ c $ in the mean function, the signal variance $ {\sigma}_f^2 $ and characteristic length scales $ \mathbf{l} $ in the covariance kernel, there is the noise variance $ {\sigma}_y^2 $ that reflects the noise level $ \varepsilon $ that is introduced independent of the objective function $ f\left(\mathbf{x}\right) $ , such that $ \theta =\left\{c,{\sigma}_f^2,\mathbf{l},{\sigma}_y^2\right\} $ .

The Gaussian process can be used to make predictions along with corresponding uncertainty quantification by computing the posterior predictive distribution (7). For $ {n}_{\ast } $ test points $ {\mathbf{X}}_{\ast } $ , it can be computed as the multivariate normal distribution conditional on the training data $ {\mathcal{D}}_n $

(7) $$ f\left({\mathbf{X}}_{\ast}\right)\hskip0.3em \mid \hskip0.3em {\mathcal{D}}_n,{\mathbf{X}}_{\ast}\sim \mathcal{N}\left({\mu}_n\left({\mathbf{X}}_{\ast}\right),{\sigma}_n^2\left({\mathbf{X}}_{\ast}\right)\right) $$

(8) $$ {\mu}_n\left({\mathbf{X}}_{\ast}\right)=K\left({\mathbf{X}}_{\ast },{\mathbf{X}}_n\right){\left[K\left({\mathbf{X}}_n,{\mathbf{X}}_n\right)+{\sigma}_y^2I\right]}^{-1}\left(\mathbf{y}-m\left({\mathbf{X}}_n\right)\right)+m\left({\mathbf{X}}_{\ast}\right) $$

(9) $$ {\sigma}_n^2\left({\mathbf{X}}_{\ast}\right)=K\left({\mathbf{X}}_{\ast },{\mathbf{X}}_{\ast}\right)-K\left({\mathbf{X}}_{\ast },{\mathbf{X}}_n\right){\left[K\left({\mathbf{X}}_n,{\mathbf{X}}_n\right)+{\sigma}_y^2I\right]}^{-1}K\left({\mathbf{X}}_n,{\mathbf{X}}_{\ast}\right), $$

where $ m\left({\mathbf{X}}_{\ast}\right) $ is the mean vector of length $ {n}_{\ast } $ over all test inputs and $ K\left({\mathbf{X}}_{\ast },{\mathbf{X}}_n\right) $ , $ K\left({\mathbf{X}}_n,{\mathbf{X}}_{\ast}\right) $ , and $ K\left({\mathbf{X}}_{\ast },{\mathbf{X}}_{\ast}\right) $ are the covariance matrices of sizes $ {n}_{\ast}\times n $ , $ n\times {n}_{\ast } $ , and $ {n}_{\ast}\times {n}_{\ast } $ between the training inputs $ {\mathbf{X}}_n $ and the test inputs $ {\mathbf{X}}_{\ast } $ , respectively.

3.1.2. Acquisition functions

Acquisition functions guide the sequential selection of candidate points by quantifying if a particular input point is likely to be a good new candidate point and thus should be evaluated by the objective function $ f\left(\mathbf{x}\right) $ . This is achieved by computing and maximizing an acquisition criterion $ \alpha \left(\cdot \right) $ based on the posterior distribution of the Gaussian process and the available training data $ {\mathcal{D}}_n $ . The criterion’s exact form depends on the individual acquisition function; however, most acquisition functions have one property in common—the exploration-exploitation trade-off (Shahriari et al. Reference Shahriari, Swersky, Wang, Adams and De Freitas2015; Frazier Reference Frazier2018). Exploration can be defined as choosing candidate points from areas with high uncertainty, where no training data points were observed, and thus, little information is available. Conversely, exploitation is defined as selecting candidate points from areas with a high predictive mean, that is, points that are close to high training data points. To understand the importance of balancing exploration and exploitation, consider the extreme cases of pure exploration and pure exploitation. For the former, only points with the highest uncertainty would be selected. While this minimizes the uncertainty of the Gaussian process, it is not a sample-efficient approach as information about high predictive means is disregarded, and areas with the best-observed outputs are avoided. For the latter, the algorithm will blindly follow the best-performing points and never explore other areas. The algorithm will probably converge toward the first optimum it discovers, making it prone to getting stuck in a local optimum. Thus, a hybrid solution that combines exploration and exploitation is beneficial.

While related work (Krause and Ong Reference Krause and Ong2011) focuses on a theoretical investigation of UCB (Srinivas et al. Reference Srinivas, Krause, Kakade and Seegre2010), this study empirically investigates EI (Jones et al. Reference Jones, Schonlau and Welch1998) and shows its superior performance on two test functions. One of the main advantages of EI over UCB is that it does not require the choice of a value for the trade-off parameter that balances exploration and exploitation. Finding optimal values for this parameter that perform well for specific problems is not trivial, although theoretical guarantees have been derived (Srinivas et al. Reference Srinivas, Krause, Kakade and Seegre2010). EI has been investigated empirically (Jones et al. Reference Jones, Schonlau and Welch1998; Jones Reference Jones2001; Brochu et al. Reference Brochu, Cora and de Freitas2010), its convergence has been studied (Bull Reference Bull2011), and it has shown efficiently in practice (Snoek et al. Reference Snoek, Larochelle and Adams2012). Although improvement-based methods are among the most popular and well-studied acquisition functions (Ament et al. Reference Ament, Daulton, Eriksson, Balandat and Bakshy2023), other acquisition functions, such as knowledge gradient (Frazier et al. Reference Frazier, Powell and Dayanik2009; Frazier Reference Frazier2018) and information-based approaches, such as entropy search (Villemonteix et al. Reference Villemonteix, Vázquez and Walter2006; Hennig and Schuler Reference Hennig and Schuler2011), predictive entropy search (Hernández-Lobato et al. Reference Hernández-Lobato, Gelbart, Hoffman, Adams and Ghahramani2015), and max-value entropy search (Wang and Jegelka Reference Wang and Jegelka2017), are also available. While we focus on EI and UCB here, our method could be used with any choice of acquisition function.

EI (Jones et al. Reference Jones, Schonlau and Welch1998) is an improvement-based acquisition function that aims to find candidates that perform better than a defined target, typically the best available training data point. EI is defined as

(10) $$ {\alpha}_{\mathrm{EI}}\left({\mathbf{X}}_{\ast}\right)=\left\{\begin{array}{ll}\left({\mu}_n\left({\mathbf{X}}_{\ast}\right)-{y}^{best}\right)\Phi (z)+{\sigma}_n\left({\mathbf{X}}_{\ast}\right)\phi (z)& \mathrm{if}\;{\sigma}_n\left({\mathbf{X}}_{\ast}\right)>0\\ {}0& \mathrm{if}\;{\sigma}_n\left({\mathbf{X}}_{\ast}\right)=0\end{array}\right., $$

where $ z=\frac{\mu_n\left({\mathbf{X}}_{\ast}\right)-{y}^{best}}{\sigma_n\left({\mathbf{X}}_{\ast}\right)} $ , $ {\mu}_n\left(\cdot \right) $ and $ {\sigma}_n\left(\cdot \right) $ are the mean and the standard deviation of the Gaussian processes’ posterior predictive distribution (7), $ {y}^{best} $ is the current best observation, and $ \Phi \left(\cdot \right) $ and $ \phi \left(\cdot \right) $ are the cumulative distribution function and probability density function of the standard normal distribution $ \mathcal{N}\left(0,1\right) $ .

Although EI can never be nonpositive mathematically, it can become zero numerically when computed due to the floating-point precision of the programming language. This results in flat areas where the EI is zero and cannot be optimized correctly. To prevent numerically vanishing values, LogEI was proposed by Ament et al. (Reference Ament, Daulton, Eriksson, Balandat and Bakshy2023) as

(11) $$ {\alpha}_{\mathrm{LogEI}}\left({\mathbf{X}}_{\ast}\right)=\left\{\begin{array}{ll}{\log}_h\left(\frac{\mu_n\left({\mathbf{X}}_{\ast}\right)-{y}^{best}}{\sigma_n\left({\mathbf{X}}_{\ast}\right)}\right)+\log \left({\sigma}_n\left({\mathbf{X}}_{\ast}\right)\right)& \mathrm{if}\;{\sigma}_n\left({\mathbf{X}}_{\ast}\right)>0\\ {}0& \mathrm{if}\;{\sigma}_n\left({\mathbf{X}}_{\ast}\right)=0\end{array}\right., $$

where

(12) $$ {\log}_h(z)=\left\{\begin{array}{ll}\log \left(\phi (z)+z\Phi (z)\right)& \mathrm{if}\;z>-1\\ {}-{z}^2/2-{c}_1+\log \hskip-0.1em 1\mathrm{mexp}\left(\log \left(\mathrm{erfcx}\left(-z/\sqrt{2}\right)|z|\right)+{c}_2\right)& \mathrm{if}-1/\sqrt{\unicode{x025B}}<z\le -1\\ {}-{z}^2/2-{c}_1-2\log \left(|z|\right)& \mathrm{if}\;z\le -1/\sqrt{\unicode{x025B}}\end{array}\right., $$

where $ {c}_1=\log \left(2\pi \right)/2 $ , $ {c}_2=\log \left(\pi /2\right)/2 $ , $ \unicode{x025B} $ is the numerical precision, and $ \log \hskip-0.22em 1\mathrm{mexp} $ and $ \mathrm{erfcx} $ are stable implementations of $ \log \left(1-\exp (z)\right) $ and $ \exp \left({z}^2\right)\operatorname{erfc}(z) $ , respectively, where erfc is the complementary error function.

The UCB (Srinivas et al. Reference Srinivas, Krause, Kakade and Seegre2010) is an optimistic acquisition function and assumes the uncertainty of the posterior Gaussian process to be true to a predefined level. It can be computed as

(13) $$ {\alpha}_{\mathrm{UCB}}\left({\mathbf{X}}_{\ast}\right)={\mu}_n\left({\mathbf{X}}_{\ast}\right)+\sqrt{\beta }{\sigma}_n\left({\mathbf{X}}_{\ast}\right), $$

where $ \beta $ is a predefined trade-off parameter that can be set for each iteration of the Bayesian optimization algorithm, it can be kept constant for the full optimization campaign or be varied at each iteration (Srinivas et al. Reference Srinivas, Krause, Kakade and Seegre2010). Srinivas et al. (Reference Srinivas, Krause, Kakade and Seegre2010) investigated some theoretical properties of $ \beta $ , while Diessner et al. (Reference Diessner, O’Connor, Wynn, Laizet, Guan, Wilson and Whalley2022) investigated how different values for $ \beta $ affect the optimization. As the acquisition functions in this study are deterministic, they can be maximized with a deterministic optimizer, such as L-BFGS-B (Zhu et al. Reference Zhu, Byrd, Lu and Nocedal1997).

3.2. Changing environmental conditions

The Bayesian optimization algorithm in Algorithm 1 assumes that all parameters influencing the output can be controlled. However, in many cases, when optimizing physical experiments, variables will be present that influence the output but cannot be controlled. This article refers to these uncontrollable variables as environmental variables as they are externally given within the ambient environment of the experiment. Examples of such uncontrollable variables are temperature, humidity and wind speed. This section presents an extension to Algorithm 1 that allows the inclusion of uncontrollable environmental variables in the optimization process. The extension can be broken down into three parts as highlighted in Algorithm 2.

The main modification to Algorithm 1 concerns the surrogate modeling. The basic Bayesian optimization algorithm fits a surrogate model over all controllable variables. Environmental variables are not included and are assumed to be fixed over the full optimization process or are irrelevant to the output. Algorithm 2 does not make this assumption and includes all controllable parameters $ {\mathbf{x}}_C $ and environmental variables $ {\mathbf{x}}_E $ in its surrogate model. The inputs $ {\mathbf{X}}_n $ of the training data $ {\mathcal{D}}_n=\left\{{\mathbf{X}}_n,{\mathbf{y}}_n\right\} $ are extended from $ {\mathbf{X}}_n=\left\{{\mathbf{X}}_{n,C}\right\} $ to $ {\mathbf{X}}_n=\left\{{\mathbf{X}}_{n,E},{\mathbf{X}}_{n,C}\right\} $ .

The second extension regards the computation of the next candidate point, specifically, the maximization of the acquisition function. While in Algorithm 1 all parameters are assumed to be controllable and the acquisition function can be maximized over all parameters $ {\max}_{{\mathbf{x}}_C}\alpha \left({\mathbf{x}}_C\right) $ , Algorithm 2 must differentiate between the controllable parameters $ {\mathbf{x}}_C $ and environmental variables $ {\mathbf{x}}_E $ . The uncontrollable variables are given by the environment and can only be measured but not manipulated. Hence, the maximization of the acquisition function is broken down into two steps. First, the environmental variables are measured. This gives values for the uncontrollable inputs for the next candidate point $ {\mathbf{x}}_{n+1,E} $ . Second, the acquisition function is maximized conditional on these values for the environmental inputs $ {\max}_{{\mathbf{x}}_C\mid {\mathbf{x}}_E}\alpha \left({\mathbf{x}}_C\right) $ resulting in the controllable inputs for the next candidate point $ {\mathbf{x}}_{n+1,C} $ . Conditional maximization essentially means that the environmental variables $ {\mathbf{x}}_E $ are treated as fixed for the maximization of the acquisition function for one iteration. This assumes that the environmental variables do not change significantly from the time of measuring until the evaluation of the new candidate point $ {\mathbf{x}}_{n+1} $ . This assumption should be realistic for most experiments, as one iteration of the Bayesian optimization loop—that is, measuring the environmental variables, fitting a Gaussian process and optimizing the acquisition function conditional on the measurements—takes only a few seconds (see Diessner et al. (Reference Diessner, Wilson and Whalley2023) for runtimes of different Bayesian optimization packages). However, issues could arise when working with environmental variables that change rapidly. Thus, Section 5 investigates the influence of different fluctuation rates. The new candidate point $ {\mathbf{x}}_{n+1} $ is then defined as a combination of the measurements for the environmental variables $ {\mathbf{x}}_{n+1,E} $ and the results of the maximization of the acquisition function $ {\mathbf{x}}_{n+1,C} $ .

The last adjustment to Algorithm 1 focuses on generating the training data. Usually, training data are produced using a space-filling design, such as a Latin hypercube (McKay et al. Reference McKay, Beckman and Conover1979). Under the assumption that all parameters can be controlled, the experiment can be repeated for each training data point to observe its output. However, the modified Bayesian optimization algorithm includes uncontrollable variables in its computation. Thus, it is impossible to evaluate any arbitrary combination of inputs as it is limited by the current measurement of the environmental variables. To resolve this issue, Algorithm 2 uses one training data point $ {\mathbf{x}}_0 $ instead of multiple points generated from a space-filling design. The initial training data are restricted to a single point, while the second data point is already computed via Bayesian optimization. The first data point is generated by taking measurements for the environmental variables $ {\mathbf{x}}_{0,E} $ and randomly selecting values for the controllable parameter $ {\mathbf{x}}_{0,C} $ . These inputs are then evaluated, resulting in a complete training inputs-output pair $ {\mathcal{D}}_0=\left\{{\mathbf{x}}_0,{y}_0\right\} $ , followed by the first Bayesian optimization loop.

Algorithm 2 Modified Bayesian optimization algorithm with environmental conditions.

In contrast to Krause and Ong (Reference Krause and Ong2011), the proposed approach does not assume different covariance structures for controllable and environmental variables. When working with experiments and simulators, the underlying objective function is generally unknown or too complex to compute directly. Even with expert knowledge, there might be insufficient information about these black boxes to confidently assume a linear or additive structure for the environmental variable (Frazier Reference Frazier2018). Hence, providing the surrogate model with enough flexibility to estimate the covariance structure is essential. This can be achieved using the Matérn kernel for controllable and environmental variables—or the radial basis function kernel for very smooth objective functions.

Figure 2 illustrates the conditional variable optimization on a two-dimensional problem with one uncontrollable variable $ {x}_1 $ and one controllable parameter $ {x}_2 $ . The upper-left plot shows the true output of the objective function, where yellow areas indicate high function values and blue areas indicate low function values. The goal is to find the optimal value for the controllable parameter (y-axis) that maximizes the output for any uncontrollable variable (x-axis) value. The upper-right plot shows the predictive mean of a Gaussian process fitted to 20 training data points (black crosses). Following Algorithm 2, a measurement (red dashed line) of the uncontrollable variable results in $ {x}_1=-0.5 $ . The next iteration of the optimization loop is performed conditional on this measurement. The lower-left plot shows the predictive mean of the Gaussian process for $ {x}_1=-0.5 $ . The conditional optimization takes a slice from the full surrogate model. It reduces the two-dimensional optimization problem to a one-dimensional problem for each optimization step, where only the controllable parameters are considered. However, the information gained from the data is shared between each iteration. Notice that no training points lie on the measurement line; however, the model uses the available training points to inform its prediction. If Algorithm 1 were used, the uncontrollable input would be assumed to be fixed for the entire optimization loop, and the optimization process would need repeating for each value of $ {x}_1 $ . The lower-right plot extends the lower-left plot by adding the uncertainty from the Gaussian process and the acquisition function. The optimal value of the controllable input $ {x}_2 $ is found by maximizing the acquisition function, and the new candidate point is a combination of this maximum and the measurement taken for the uncontrollable variable $ {x}_1 $ . The candidate point is observed and added to the training data to be used in the next iteration of the optimization loop.

Figure 2. Maximization of a two-dimensional problem with one environmental variable $ {x}_1 $ and one controllable variable $ {x}_2 $ . Yellow areas indicate high outputs, and dark blue areas indicate low outputs. Upper-left: True objective function. Upper-right: Prediction of a Gaussian process with a measurement taken for the following conditional optimization step. Lower-left: Gaussian process prediction for optimization conditional on the measurement. Lower-right: Bayesian optimization step conditional on the measurement.

4. Simulations

This section introduces two synthetic test functionsFootnote ¹—the Levy function and the Hartmann function—and applies the Bayesian optimization algorithm with environmental conditions presented in Section 3.2. Simulations for both problems are run for 100 evaluations and are repeated 30 timesFootnote ² to validate the robustness of Algorithm 2. This decreases the risk that results are influenced by the method’s inherent randomness, for example, the randomly sampled training points that initialize the algorithm. Both problems assume one uncontrollable variable whose value is provided by a random walk at each iteration. In the simulations, each step of the random walk adds a sample from a uniform distribution $ \mathcal{U} $ to the previous value of the uncontrollable variable, such that

(14) $$ {\mathbf{x}}_{n,E}={\mathbf{x}}_{n-1,E}+{\mathcal{U}}_{\left[-\mathbf{a},\mathbf{a}\right]}, $$

where $ \mathbf{a} $ is a vector of small predefined constants that provide the minimal and maximal change of the environmental variables from one iteration to the next. This uniform assumption represents the natural fluctuation of the uncontrollable variables encountered in a real-world application, for example, changes in temperature, humidity, and wind speed. It further allows the investigation of uncontrollable variables with different fluctuation levels by increasing or decreasing the constants in $ \mathbf{a} $ as discussed in Section 5. Assuming another distribution for the constants in $ \mathbf{a} $ , such as a Normal distribution, is an alternative to this approach.

Subsequent sections analyze the performance of Algorithm 2 by comparing predictions from the Gaussian process models $ {\mu}_n\left(\mathbf{x}\right) $ to the actual optimal values $ f\left(\mathbf{x}\right) $ . In both cases, results are obtained by maximizing the Gaussian process prediction and the true objective function conditional on identical test values $ {\mathbf{x}}_E^{\prime } $ for the uncontrollable variable. These test values are sampled from a maximin Latin hypercube design (McKay et al. Reference McKay, Beckman and Conover1979) within the observed domain of the uncontrollable variable $ \left[\min \left({\mathbf{x}}_{N,E}\right),\max \left({\mathbf{x}}_{N,E}\right)\right] $ considered by the algorithm. We also call this observed domain the effective domain. This ensures a fair comparison by avoiding predictions outside the effective domain requiring extrapolation. Extrapolation with Gaussian processes generally means that predictions default to the prior mean function. In cases where extrapolation cannot be avoided, making the prior mean function as informative as possible—for example, by going beyond zero and constant mean functions with polynomial and trigonometric mean functions—can improve results significantly (Planas et al. Reference Planas, Oune and Bostanabad2020). To score the performance of the algorithm, the mean absolute percentage error $ \mathrm{MAPE}\left({\mu}_n\left({\mathbf{X}}_i^{\prime}\right),f\left({\mathbf{X}}_i^{\prime}\right)\right)=\frac{1}{m}{\sum}_{i=1}^m\mid \frac{\Big({\mu}_n\left({\mathbf{x}}_i^{\prime}\right)-f\left({\mathbf{x}}_i^{\prime}\right)}{f\left({\mathbf{x}}_i^{\prime}\right)}\mid $ between the Gaussian process’ predictions and the truths are computed for all $ m $ test points $ {\mathbf{X}}^{\prime } $ .

The acquisition criterion conditional on values of the uncontrollable variables $ {\mathbf{x}}_{n,E} $ is maximized with the SLSQP algorithm (Kraft Reference Kraft1994) using multiple starts. For this strategy, 100 points are sampled from a maximin Latin hypercube design (McKay et al. Reference McKay, Beckman and Conover1979) and evaluated by the acquisition criterion. The best 20 points are then used to initialize the SLSQP algorithm, and only the best result is used as the solution for the optimization problem. The multiple starts aim to reduce the risk of converging toward a local optimum instead of the desired global optimum of the acquisition function.

4.1. The two-dimensional Levy function

The Levy function

$$ f\left(\mathbf{x}\right)={\sin}^2\left(\pi {w}_1\right)+{\left({w}_1-1\right)}^2\left[1+10{\sin}^2\left(\pi {w}_1+1\right)\right]+{\left({w}_2-1\right)}^2\left[1+{\sin}^2\left(2\pi {w}_2\right)\right], $$

where $ {w}_i=1+\frac{x_i-1}{4} $ , for $ i=1,2 $ , is a two-dimensional function with two input parameters $ {\mathbf{x}}_1 $ and $ {\mathbf{x}}_2 $ . Although often used as a minimization problem to find the global minimum $ f\left({\mathbf{x}}^{\ast}\right)=0 $ at $ {\mathbf{x}}^{\ast }=\left(1,1\right) $ for the input space $ {\left[-10,10\right]}^2 $ , we choose the bounds of the parameters as $ \left[-\mathrm{7.5,7.5}\right] $ and $ \left[-10,10\right] $ , respectively to create a maximization problem that is better suited for testing Algorithm 2. The controllable parameter is restricted from its usual range of $ \left[-10,10\right] $ to $ \left[-\mathrm{7.5,7.5}\right] $ to avoid the steep ridges of the Levy function at $ {x}_1=-10 $ and $ {x}_1=10 $ that push the optima toward the outer bounds of the parameter space. The function given in Figure 3 shows a clear ridge at values of about $ {x}_1=-6 $ for the controllable parameter $ {x}_2 $ for all values of the uncontrollable variable $ {x}_1 $ . The simulations use $ a=1.5 $ as the uniform distribution constant of the random walk for the uncontrollable variable $ {x}_2 $ . For this relatively simple two-dimensional function, the constant $ a $ was chosen to change by no more than 7.5% of the uncontrollable parameter’s range in either direction per iteration. This assumes that the environmental conditions do not change too much between iterations. However, the effect of setting $ a $ to different values is investigated and discussed in Section 5.

Figure 3. Two-dimensional negated Levy function with one controllable parameter $ {x}_1 $ bounded by $ \left[-\mathrm{7.5,7.5}\right] $ and one uncontrollable variable $ {x}_2 $ bounded by $ \left[-10,10\right] $ .

The upper-left plot in Figure 4 shows the performance of Algorithm 2 as the mean absolute percentage error between the maximum of the Gaussian process prediction and the maximum of the true objective function conditional on 25 test values of the uncontrollable variable. Three alternatives for the acquisition function—EI, LogEI, and UCB introduced in Section 3.1.2—are compared to a benchmark where values for the controllable variable were selected randomly. Three different values (4, 8, 16) for the trade-off parameter $ \beta $ of UCB are considered in Figure 5. While results for all parameter values are comparable, UCB with trade-off parameter $ \beta =8 $ performs slightly better in terms of mean performance and variance between runs. Thus, comparisons with EI and LogEI use $ \beta =8 $ . The solid lines indicate the mean performance, while the shaded areas indicate the 95% confidence interval over the 30 replications. Each replication’s mean absolute percentage error is computed for every 10 evaluations for the same 25 test points $ {\mathbf{X}}_E^{\prime } $ of the environmental variable. The test values are sampled from a Latin hypercube bounded by the minimal and maximal value of the uncontrollable variable after 100 function evaluations. The mean absolute percentage error starts just below 0.8 for all alternatives. While the improvement-based algorithms (EI and LogEI) performed better than the random benchmark, the algorithm using UCB performs worse. After 100 function evaluations, EI and LogEI have a mean absolute percentage error of 0.08 and 0.06, respectively, and improve upon the random benchmark (0.17). UCB only achieves a mean absolute percentage error of 0.24. Moreover, the optimistic strategy shows large confidence intervals, indicating that the method is not robust. Altering the trade-off parameter $ \beta $ did not improve this result, as illustrated in the left plot in Figure 5.

Figure 4. Upper row: Means (lines) and $ 95 $ % confidence intervals (shaded areas) of the mean absolute percentage error between Gaussian process prediction and truth over $ 30 $ replications. Lower row: Difference between algorithms and random benchmark after $ 100 $ function evaluations for each of the $ 30 $ replications. Two-dimensional Levy function with one uncontrollable parameter on the left and six-dimensional Hartmann function with one uncontrollable parameter on the right.

Figure 5. Comparison of different trade-off parameters $ \beta $ for the UCB acquisition function. Means (lines) and 95% confidence intervals (shaded areas) of the mean absolute percentage error between Gaussian process prediction and truth over 30 replications. Two-dimensional Levy function with one uncontrollable parameter on the left and six-dimensional Hartmann function with one uncontrollable parameter on the right.

The lower-left plot presents the difference between the mean absolute percentage error of the random benchmark and the three alternative acquisition functions after 100 evaluations. This difference is computed for each of the 30 replications, and distributions of the differences are plotted for the three different acquisition functions. Negative values indicate replications where the benchmark resulted in superior solutions, while positive values indicate that the given version of Algorithm 2 performed better than the benchmark. Although no alternative is better than the benchmark for every single replication of the 30 total replications, there is a clear difference between the improvement-based and the optimistic acquisition functions. Indeed, the mean of UCB (displayed by the horizontal line toward the center of each violin plot) is the only one worse than zero. This means that the random benchmark outperforms UCB on average. Additionally, the plot mirrors the lack of robustness discovered earlier by the large spread in differences. While the method performs better for some replications than the random benchmark, it performs much worse for others. A Mann–Whitney U test was performed to determine what alternatives perform differently from the random benchmark to a 1% significance level. The test returned a $ P $ -value of 0.0 for the improvement-based algorithms and 0.37 for the UCB. This rejects the hypothesis that methods perform equally well for the improvement-based methods, reinforcing the results from the visual analysis that EI and LogEI perform significantly better than the random benchmark. This cannot be said for the optimistic method—the test cannot reject the null hypothesis, indicating that UCB does not perform significantly differently from a random approach.

4.2. The six-dimensional Hartmann function

The negatedFootnote ³ Hartmann function

$$ f\left(\mathbf{x}\right)=\sum \limits_{i=1}^4{\alpha}_i\exp \left(-\sum \limits_{j=1}^6{\mathbf{A}}_{ij}{\left({x}_j-{\mathbf{P}}_{ij}\right)}^2\right), $$

where

$$ {\displaystyle \begin{array}{l}\alpha ={\left(\mathrm{1.0,1.2,3.0,3.2}\right)}^T,\\ {}\boldsymbol{A}=\left(\begin{array}{cccccc}10.00& 3.00& 17.00& 3.50& 1.70& 8.00\\ {}0.05& 10.00& 17.00& 0.10& 8.00& 14.00\\ {}3.00& 3.50& 1.70& 10.00& 17.00& 8.00\\ {}17.00& 8.00& 0.05& 10.00& 0.10& 14.00\end{array}\right),\mathrm{and}\\ {}\boldsymbol{P}={10}^{-4}\left(\begin{array}{cccccc}1312& 1696& 5569& 124& 8283& 5886\\ {}2329& 4135& 8307& 3736& 1004& 9991\\ {}2348& 1451& 3522& 2883& 3047& 6650\\ {}4047& 8828& 8732& 5743& 1091& 381\end{array}\right),\end{array}} $$

is a six-dimensional function with six input parameters $ {x}_1 $ , $ {x}_2 $ , $ {x}_3 $ , $ {x}_4 $ , $ {x}_5 $ , and $ {x}_6 $ that are evaluated on the hypercube $ {\left(0,1\right)}^6 $ . It has six local maxima and one global maximum with $ f\left({\mathbf{x}}^{\ast}\right)=3.32 $ at $ {\mathbf{x}}^{\ast }=\left(\mathrm{0.20,0.15,0.48,0.28,0.31,0.66}\right) $ . The simulations use $ a=0.05 $ as the uniform distribution constant of the random walk for the environmental variable $ {x}_6 $ . As this six-dimensional function is more complex than the Levy function, the constant $ a $ was chosen to change by a maximum of 5% of the uncontrollable parameter’s range in either direction, that is, the environmental conditions do not change too much between iterations. Section 5 sets $ a $ to 10%, 25%, 50%, and 100% and investigates their effect.

The upper-right plot in Figure 4 shows the performance of Algorithm 2 over the 30 repeats for the same three acquisition function as for the Levy function and compares it to the random benchmark. The mean absolute percentage error starts between 0.85 and 0.90 for all four algorithms. After 100 function evaluations, the mean absolute percentage error between the prediction of the Gaussian process and the true objective value for Algorithm 2 using EI and UCB with $ \beta =8 $ are 0.07 and 0.06, respectively—much better than the random benchmark with 0.24. However, the algorithm using LogEI with a mean absolute percentage error of 0.18 after 100 evaluations performs only slightly better than the benchmark and significantly worse than the other versions of Algorithm 2. The large 95% confidence intervals for LogEI indicate that Algorithm 2 with LogEI is not robust in this case, making it less reliable than EI and UCB.

The violin plots of the difference between the mean absolute percentage error of the benchmark and the three variations of Algorithm 2 for all 30 replications on the lower-right in Figure 4 indicate that EI performs the best with almost all replications better than the random benchmark. LogEI, on the other hand, performs comparably to the benchmark on average but has a large spread, with some replication performing much worse. UCB performs similarly to EI but has an outlier that performs much worse than the random benchmark. Overall, EI is the most robust method and is not prone to outliers. Despite these differences between algorithms, results after 100 evaluations for all three algorithms are significantly different from the random results to a 1% significance level: the $ P $ -values from Mann–Whitney U tests are 0.0 for the EI and the UCB version of Algorithm 2 and 0.007 LogEI. While this shows that the results of all methods are significantly different from the benchmark, only EI and UCB perform better than the benchmark indicated by the better average mean absolute percentage error.

Figure 5 shows the results for different values of the trade-off parameter $ \beta $ . While the average performance is very similar after 100 evaluations, there is a difference in the 95% confidence intervals for the Hartmann function. However, no clear correlation is noticeable as $ \beta =8 $ performs better than $ \beta =4 $ and $ \beta =16 $ , indicating no clear trend.

5. Empirical analysis of properties

Based on the investigations of Section 4, we define ENVBO as a version of Algorithm 2 that uses the EI acquisition function. This section explores the effect of changes to five aspects of the underlying objective function or environmental conditions on ENVBO using the negated Hartmann function. The effect of adding different levels of random Gaussian noise to the function’s output, the influence of more than one uncontrollable variable, the effect of the fluctuation level (i.e., the step size $ a $ of the random walk), the impact of the variability in the uncontrollable variable, and the relationship between the algorithm performance and the effective domain size of the uncontrollable variables are investigated. An off-the-shelf version of ENVBO is available via the open-source Python package NUBO (Diessner et al. Reference Diessner, Wilson and Whalley2023).

To make the comparison fair, the 30 runs used the same 30 initial starting points and random walks. For the random walks, this is achieved by sampling changes in percentages and scaling them by the maximal step size $ a $ rather than sampling absolute values directly.

Most physical experiments in engineering can only be conducted by introducing some noise, such as measurement uncertainty, that cannot be eliminated entirely. This section replicates this situation by adding randomly generated noise $ \unicode{x025B} $ to the Hartmann function. The noisy function can be defined as $ g\left(\mathbf{x}\right)=f\left(\mathbf{x}\right)+\unicode{x025B} $ , where $ f\left(\mathbf{x}\right) $ is a deterministic negated Hartmann function from Section 4.2. The noise is sampled from a Normal distribution centered around zero with a small standard deviation $ \sigma $ , such that $ \varepsilon \sim \mathcal{N}\left(0,{\sigma}^2\right) $ . The simulations explore noise levels with $ \sigma =0.00 $ , $ \sigma =0.025 $ , $ \sigma =0.050 $ , and $ \sigma =0.100 $ . Considering the range of the Hartmann function, this corresponds to standard deviations of 0.75%, 1.5%, and 3.0% of the full output range, respectively. This means that for any of these three cases, 68.3% of the added noise values will fall between $ \pm 1\sigma $ , 95.5% fall between $ \pm 2\sigma $ and 99.7% fall between $ \pm 3\sigma $ . For $ \sigma =0.100 $ , this translates into noise values that decrease or increase the real output by up to 3.0% of the output range 68.3% of the time, by 6.0% of the output range 95.5% of the time, and by 9.0% of the output range 99.7% of the time. The plot in the upper-left of Figure 6 shows the performance of ENVBO for each of these four noise levels. The results indicate no significant difference in the average performance or the 95% confidence intervals between the four cases. Overall, the proposed method is not sensitive to adding modest noise levels.

Figure 6. Means (lines) and 95% confidence intervals (shaded areas) of the mean absolute percentage error between the predictive mean of the Gaussian process and the truth over 30 replications for the six-dimensional Hartmann function. Upper-left: Comparison of randomly added noise levels, $ \mathcal{N}\left(0,{\sigma}^2\right) $ . Upper-right: Comparison of different numbers of uncontrollable parameters $ {n}_E $ . Lower-left: Comparison of five different step sizes $ a $ for the random walk $ {\mathcal{U}}_{\left[-a,a\right]} $ added to the previous uncontrollable value. Lower-right: Comparison of uncontrollable variables with different parameter variability.

For some experiments, there might be more than one influential uncontrollable variable. The upper-right plot of Figure 6 provides results about the performance of ENVBO with one, two and three uncontrollable variables. At the same time, the overall dimensionality of the problem stays the same at $ n=6 $ . For $ {n}_E=1 $ , input six of the Hartmann function from Equation (4.2) is assumed uncontrollable, while input one is added to the uncontrollable variables for $ {n}_E=2 $ , and input one and four are added for $ {n}_E=3 $ . The number of test points used to evaluate the final Gaussian process model is increased from 25 to 50 for $ {n}_E=2 $ and 75 for $ {n}_E=3 $ . The results show that the mean absolute percentage error increases with increasing numbers of environmental variables. Particularly, the 95% confidence intervals widen significantly. This result is expected as the input space of the environmental variables grows exponentially with $ {n}_E $ and requires exponentially more training points to cover the input space equally well as lower $ {n}_E $ . Thus, more evaluations are required to achieve similar results.

Uncontrollable variables will fluctuate to different extents from one evaluation to the next, for example, when measured in physical experiments. Higher fluctuations cause big jumps in the uncontrollable variable values, while uncontrollable variables will be more stable for lower fluctuations. The lower-left plot of Figure 6 shows the results for five different fluctuation levels implemented by varying parameter $ a $ of the uniform distribution in Equation (14)—the higher $ a $ , the higher the fluctuation of the uncontrollable variable. The results show that Algorithm 2 performs better for lower $ a $ , and the performance of the final Gaussian process models decreases with increasing fluctuation.

Parameter variability indicates how much the output changes when moving along a parameter’s input axis. Uncontrollable variables with a low variability will only change the output slightly, while variables with a high variability will change the output considerably. The effect of value changes of two variables could differ considerably if they have different levels of parameter variability. As a proxy for the parameter variability, the length scales of a well-fitting Gaussian process over the full domain are considered. The length scales quantify how long the outputs are correlated when moving along the input axes (Gramacy Reference Gramacy2020). Consider, for example, a parameter with a large length scale. When this parameter value is changed by a certain amount, a relatively small change is expected in the output—provided everything else stays the same. The change in the output is expected to be more significant for a parameter with a small length scale. A Gaussian process is fitted to 2000 data points sampled from a Latin hypercube (McKay et al. Reference McKay, Beckman and Conover1979) to achieve a well-fitting model. The resulting length scales for all six parameters in order are 1.30, 1.90, 4.81, 1.48, 1.51, and 1.47. The first input has the lowest length scale, suggesting that outputs are only correlated for a short distance when moving along its axis. This means that the parameter variability of the first input is high. In contrast, the third input has the highest length scale, indicating a low parameter variability. Moving along its axis, less change is expected for the third input than for the first input. The lower-right plot of Figure 6 compares the first (high), third (low), and sixth input (medium) when chosen as the uncontrollable variable. Differences in the low and medium parameter variability are minimal. At the same time, there is some difference compared to the parameter with a high variability: the confidence interval is noticeably more significant, and the average of the mean absolute percentage error is slightly worse, especially for evaluations 30–80.

Finally, the relationship between the size of the effective domain, that is, the actual searched input space of the environmental variables, and the performance of the Gaussian process is investigated. Figure 7 uses the same data as Figure 6 but plots the effective parameter domain against the mean absolute percentage error, which provides a proxy for the performance of the Gaussian process. Each point reflects one individual replication of the 30 performed replications. The trend lines in each plot show a positive relationship between the effective domain and the mean absolute percentage error. Small effective domains generally correspond to small mean absolute percentage errors, while large effective domains generally correspond to large mean absolute percentage errors. Only the upper-right plot that plots the effective domain for different numbers of uncontrollable parameters against the mean absolute percentage error shows almost no relationship.

Figure 7. Relationships between the actual effective domain of the uncontrollable variables and the mean absolute percentage error of an individual run for the six-dimensional Hartmann function. Upper-left: Comparison of randomly added noise levels, $ \mathcal{N}\left(0,{\sigma}^2\right) $ . Upper-right: Comparison of different numbers of uncontrollable parameters $ {n}_E $ . Lower-left: Comparison of five different step sizes $ a $ for the random walk $ {\mathcal{U}}_{\left[-a,a\right]} $ added to the previous uncontrollable value. Lower-right: Comparison of uncontrollable variables with different parameter variability.

6. Application to a wind farm simulator

The power generation of wind farms is highly dependent on the wind speed, the wind direction and the placement of the individual wind turbines within a particular area or site (Grady et al. Reference Grady, Hussaini and Abdullah2005; Qureshi and Warudkar Reference Qureshi and Warudkar2023). Most of the time, interest lies in finding the optimal wind turbine positions that maximize the energy production conditional on either constant or variable wind speeds and directions (Mosetti et al. Reference Mosetti, Poloni and Diviacco1994; Chen et al. Reference Chen, Li, Jin and Song2013; Parada et al. Reference Parada, Herrera, Flores and Parada2017). Optimization algorithms such as TOPFARMFootnote ⁴ can be used when function evaluations are cheap, while methods discussed in Section 2 (Williams et al. Reference Williams, Santner and Notz2000; Swersky et al. Reference Swersky, Snoek and Adams2013; Toscano-Palmerin and Frazier Reference Toscano-Palmerin and Frazier2018) present a cost-effective alternative when function evaluations are expensive. However, these methods cannot find multiple solutions for different environmental conditions within one optimization run, which is the objective of this section. Specifically, this application aims to find the optimal wind turbine positions conditional on randomly changing wind directions, resulting in one solution for each possible wind direction. While the wind speed is assumed to be fixed in this application, it could also be randomly changing. The result would be solutions for all combinations of wind direction and wind speed.

The first row of plots of Figure 8 shows the effect of the wind direction on the local wind speed for a complex underlying terrain while the global wind speed is fixed at 6 m/s. The locations of the high local wind speeds—necessary for high-energy production—shift significantly between a wind direction of 0 and 120 degrees. While there is a band of high local wind speeds down the center of the X location for the latter, it roughly rotates 90 degrees for the former. The plots show that ideal wind turbine positions will likely differ for the two wind directions. Another critical factor for maximizing energy production is taking the wake of the wind turbines into account (lower row of Figure 8). While the wake of the wind turbines located upstream does not affect wind turbines located downstream for a wind direction of 120, the wake of wind turbine 3 for a wind direction of 0 degrees heavily affects wind turbine 1. This shows that although a wind farm position can be ideal for one wind speed, it can be suboptimal for another.

Figure 8. Wind farm simulator. Upper row: Local wind speed over the complex terrain for a wind direction of 0 and 120 degrees. Lower row: Wake of four wind turbines for a wind direction of 0 and 120 degrees. Wind speed is fixed at 6 m/s.

In this section, a fictitious site is considered with complex terrain as shown in Figure 8 on which four wind turbines have to be placed to maximize the annual energy production (AEP). The wind direction is an environmental variable that varies according to a random walk as defined in Equation (14) where 90 and 135 degrees are the lower and upper bounds, respectively, and wind directions can change by $ \pm 5 $ degrees from one iteration to the next. The global wind speed is fixed at 6 m/s. This results in a nine-dimensional problem with eight controllable parameters—one X and one Y location for each of the four wind turbines—and the wind direction as the only environmental variable. Additionally, the positioning of the wind turbines is constrained such that wind turbines have to be at least 160 meters apart to prevent the blades from colliding. Simulations are performed with PyWake (Pedersen et al. Reference Pedersen, Forsting, van der Laan, Riva, Romàn, Risco, Friis-Møller, Quick, Christiansen, Rodrigues, Olsen and Réthoré2023) and use Vestas V80 wind turbines that can produce 2 MW of energy and have a blade diameter of 80 meters. Simulating the power generation of wind farms is a complex task that requires an expensive simulator that accounts for many different parameters. However, the main objective of this example is to illustrate the performance of ENVBO and enable other researchers to replicate and validate this work. Thus, a relatively simple simulator was chosen to make the results easily reproducible so that ENVBO can be used with confidence on more expensive physical experiments and computer simulators.

ENVBO—a version of Algorithm 2 with EI as its acquisition function—is run for a function evaluation budget of 200 and benchmarked against regular Bayesian optimization (Algorithm 1) and the SLSQP optimization algorithm (Kraft Reference Kraft1994). Regular Bayesian optimization, referred to as BO in the following paragraphs, uses EI as its acquisition function to make comparisons fair. While ENVBO can return one solution for each wind direction after one optimization run, BO and the SLSQP algorithm must be run for each possible wind direction. To compare the algorithms, four wind directions are chosen—90, 105, 120, and 135 degrees—and BO is restricted to 50 function evaluations each to reach the same function evaluation budget as ENVBO. The function evaluations of the SLSQP algorithm cannot be restricted, and the algorithm is therefore run until convergence. BO and SLSQP are run for fixed wind directions, so they do not use the random walk in contrast to ENVBO. This makes converging toward a solution easier, as the algorithms can control all influential variables. ENVBO and BO were implemented via the open-source package NUBO (Diessner et al. Reference Diessner, Wilson and Whalley2023), and an off-the-shelf version of ENVBO is available at www.nubopy.com. The proposed SLSQP algorithm was implemented via the SciPy package (Virtanen et al. Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett, Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson, Carey, Polat, Feng, Moore, VanderPlas, Laxalde, Perktold, Cimrman, Henriksen, Quintero, Harris, Archibald, Ribeiro, Pedregosa and van Mulbregt2020). Furthermore, all code for optimizing the wind farm simulator is available at https://github.com/mikediessner/environmental-conditions-BO.

Figure 9 shows the results for all three algorithms and all four wind directions. The dashed circles around the wind turbine positions indicated with crosses represent the placement constraint—no wind turbine of one color can be placed within the dashed circle of another wind turbine with the same color. For a wind direction of 90 degrees, ENVBO performs best with an annual power production of 4.20 GWh, followed by BO and SLSQP with 3.51 and 2.42 GWh, respectively. This is a 20% improvement over BO and a 74% improvement over SLSQP. SLSQP performs much worse than ENVBO and places at least one turbine in a subpar area with low local wind speeds, possibly due to converging toward a local maximum. While BO places three turbines in areas with high local wind speeds, it places one within a suboptimal area. For a wind direction of 105 degrees, ENVBO with 4.70 GWh performs significantly better than BO and SLSQP, achieving 0.91 and 1.86 GWh less power generation, respectively—a 24% and 65% improvement. A similar result is achieved for a wind direction of 120 degrees. ENVBO outperforms BO and SLSQP by 0.81 and 2.4 GWh or 19% and 88%, respectively. The results of the three strategies are closest for a wind direction of 135 degrees. This is the only instance where a benchmark beats ENVBO—in this case, BO with an AEP of 2.88 GWh. ENVBO achieves 0.29 GWh (10%) less AEP but still outperforms SLSQP by 0.13 GWh (5%). However, the differences are much smaller than for any of the first three wind directions. Considering Figure 9, ENVBO is the only strategy that consistently places all four wind turbines in the band of high local wind speeds located down the center of the X location. At least one turbine is placed off to one side of this band for all other methods.

Figure 9. Annual energy production and placement of four wind turbines with spacing constraints. ENVBO (Algorithm 2) is benchmarked against SLSQP and BO (Algorithm 1).

While the results in Figure 9 show that ENVBO outperforms BO and SLSQP in almost all instances, they do not consider the different numbers of function evaluations required by each algorithm. Figure 10 plots the AEP and the number of function evaluations against the wind direction. The lower plot shows that SLSQP uses the most function evaluations by far, with 323, 294, 324, and 301 evaluations for the four wind directions—a total budget of 1152 function evaluations. In contrast, ENVBO and BO both use 200 function evaluations in total. BO divides this budget equally over the four wind directions, allocating 50 evaluations per wind direction, while they are distributed via a random walk for ENVBO. Each of ENVBO’s bins of the lower plot only contains two to 24 function evaluations. Compared to SLSQP, both Bayesian optimization algorithms are much more sample-efficient and use less than 20% of its evaluation budget. ENVBO uses less than half the evaluation for the four wind directions compared to BO, but it outperforms BO for three of the four wind directions, as shown in the upper plot. This shows the advantage of a global surrogate model that is fit to all controllable and environmental variables. The Gaussian process uses all available information and can model the effect of the environmental conditions.

Figure 10. Annual energy production and number of function evaluations conditional on the wind direction. ENVBO (Algorithm 2) is benchmarked against SLSQP and BO (Algorithm 1).

Although the improved performance of ENVBO compared to BO and SLSQP is already valuable, the main advantage lies in ENVBO’s capability to predict a solution for each possible wind direction. This is illustrated by the 51 solutions given for ENVBO in the upper plot of Figure 10. SLSQP and BO can only give solutions for the specific wind direction for which they were run. To achieve similar results, SLSQP and BO need to be rerun for each wind direction, which would multiply the required function evaluations many times. ENVBO uses the available function evaluation budget much more effectively and is a more sample-efficient and cost-effective approach.

Overall, ENVBO finds solutions over the whole range of the environmental variable, while in most cases performing up to 88% better than algorithms that focus on one fixed environmental variable at a time. This is even true for ranges of the environmental variable that are only explored briefly as the algorithm learns the effect of the environmental variable from adjacent areas. Furthermore, ENVBO uses only a small fraction of the evaluation budget of the two benchmarks, making it particularly beneficial for optimizing computer simulators and physical experiments that can be very expensive to run in engineering.