# Constrained Bayesian optimization with merit functions

Jingyi Wang
Lawrence Livermore National Laboratory
Livermore, CA 94550
<EMAIL_ADDRESS>
Cosmin G. Petra
Lawrence Livermore National Laboratory
Livermore, CA 94550
<EMAIL_ADDRESS>
&J. Luc Peterson
Lawrence Livermore National Laboratory
Livermore, CA 94550
<EMAIL_ADDRESS>

###### Abstract

Bayesian optimization is a powerful optimization tool for problems where
native first-order derivatives are unavailable. Recently, constrained Bayesian
optimization (CBO) has been applied to many engineering applications where
constraints are essential. However, several obstacles remain with current CBO
algorithms that could prevent a wider adoption. We propose CBO algorithms
using merit functions, such as the penalty merit function, in acquisition
functions, inspired by nonlinear optimization methods, e.g., sequential
quadratic programming. Merit functions measure the potential progress of both
the objective and constraint functions, thus increasing algorithmic efficiency
and allowing infeasible initial samples. The acquisition functions with merit
functions are relaxed to have closed forms, making its implementation readily
available wherever Bayesian optimization is. We further propose a unified CBO
algorithm that can be seen as extension to the popular expected constrained
improvement (ECI) approach. We demonstrate the effectiveness and efficiency of
the proposed algorithms through numerical experiments on synthetic problems
and a practical data-driven engineering design problem in the field of plasma
physics.

##  1 Introduction

Bayesian optimization is a black-box optimization method that does not require
derivative information of the objective Wang et al. (2023). Instead, a
surrogate model, often a Gaussian process (GP), is constructed to replace the
objective function based on iterative sampling of the data. It has been widely
adopted in many applications such as inverse problems Wang and Zabaras (2004),
structural design Mathern et al. (2021), additive manufacturing process
simulation Wang and Papadopoulos (2023), robotics Calandra et al. (2016),
inertial confinement fusion (ICF) design Wang et al. (2024), etc. The method
has been successfully extended to multifidelity Bayesian optimization with
multifidelity surrogate models Zuluaga et al. (2013) and Gaussian processes
with independent constraints Bernardo et al. (2011).

Constrained Bayesian optimization (CBO) has attracted huge interest recently
since constraints are critical to many engineering applications. The
mathematical formulation of an inequality constrained optimization problem is

$\displaystyle\underset{\begin{subarray}{c}x\end{subarray}\in
C}{\text{minimize}}$ $\displaystyle f(x),$ (1) subject to $\displaystyle
c(x)\geq 0,$

where $C\subset\mathbb{R}^{n}$ represents the bound constraints on $x$,
$f:\mathbb{R}^{n}\to\mathbb{R}$ is the objective function, and
$c:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is the constraint functions.

One of the most successful and widely adopted CBO methods is the expected
constrained improvement (ECI) Schonlau et al. (1998); Gelbart et al. (2014);
Gardner et al. (2014), where the expected improvement of the objective is
multiplied by the probability of constraints being satisfied. An important
reason to its success is that no new statistics routines are needed in its
implementation compared to expected improvement, as the probability of
feasibility only requires the cumulative distribution function (CDF) and
probability density function (PDF) of the standard normal distribution. Its
extension to noisy setting and scalable methods can be found in Letham et al.
(2019); Eriksson and Poloczek (2021), respectively. To account for useful
information provided by infeasible samples, integrated conditional expected
improvement (IECI) Bernardo et al. (2011), and expected volume reduction (EVR)
Picheny (2014) are proposed. Other than expected improvement, information-
based acquisition functions have seen success in the constrained setting with
methods such as the predictive entropy search with constraints (PESC)
Hernández-Lobato et al. (2015). In Zhang et al. (2023), the authors propose an
active learning algorithm and a corresponding acqusition function to search
for feasible region.

Augmented Lagrangian (AL) methods are another important group of CBO methods
that draw ideas from the AL methods in nonlinear optimization Gramacy et al.
(2016). These methods construct augmented Lagrangian functions based on the
values of both the objective and the constraints. The square of a Gaussian
process leads to a new distribution that could require approximation by Monte
Carlo integration. In the latest iteration of AL methods, Slack-AL Picheny et
al. (2016) and ADMMBO .Ariafar et al. (2019), a closed form acqusition
function has been achieved, while the algorithms themselves have become more
complex. Recently, constraint violation has been proposed as an improvement on
ECI algorithms when no feasible samples are available Jiao et al. (2019). The
authors design a separate branch of the algorithm using the constraint
functions alone when no feasible points are available. In Zhao and Xu (2024),
the authors opted for a smoothed exact penalty function and choose scaled
expected improvement (SEI) as the acquisition function. However, SEI can
encounter zero regions, in which case the algorithm resorts to a mean value
reduction scheme.

Despite the recent progress, there remain several limitations for majority of
the CBO methods. First, many methods, including ECI, IECI and EVR, require
feasible starting points, which is not available for many engineering
problems. Second, some of the more recent methods that have superior reported
performances such as IECI, EVR, AL Bayesian optimization Gramacy et al. (2016)
and PESC do not have have closed forms, and therefore require computationally
expensive estimate routines in their acquisition functions. Third, the
constraints in methods such as ECI are assumed to be conditionally
independent, which can potentially hinder their abilities to find feasible
optimal solutions. Finally, some methods that enjoy closed form acqusition
functions and excellent reported performance such as ADMMBO are not readily
available in different coding languages as they require additional libraries
and more complicated implementation.

In this paper, we propose two novel merit constrained Bayesian optimization
(MCBO) algorithms using expected merit improvement (EMI) acqusition functions,
both based on the penalty merit function and expected improvement. We further
propose a unified constrained Bayesian optimization (UCBO) algorithm by
combining EMI and ECI for a unified expected constrained improvement (UECI)
acquisition function. The essence of our algorithms is the additive structure
from EMI and a relaxation from the double maximum functions on the constraint
functions to one. Our algorithms can overcome all the limitations mentioned
above. That is, it allows infeasible starting point, has closed-form
expressions, do not require constraints to be conditionally independent, and
finally, can be easily implemented wherever unconstrained Bayesian
optimization is. Moreover, the unified algorithm can also been seen as an
extension or improvement to ECI.

This paper is organized as follows. In section 2, we describe the general
Bayesian optimization method with expected improvement acquisition function
and its extension to the constrained setting via ECI. In Section 3, we
motivate and propose our MCBO and UCBO algorithms and their implementation.
Numerical experiments on three synthetic problems and a data-driven ICF design
test problem are presented in section 4 that illustrate and compare the
capabilities of the proposed methods. Conclusions are given in section 5.

##  2 Bayesian optimization

In Bayesian optimization, a subproblem of (1) is solved at each iteration,
where a surrogate function replaces $f$. A surrogate function is an
approximation of $f$ constructed based on a collection of available data
samples. It is formed and updated at each iteration to obtain a more accurate
approximation as the optimization progresses. A Gaussian process surrogate
model assumes that at sample point $x$, a multivariate jointly Gaussian
distribution describes the existing samples and the design objective $f$ as
follows:

$\displaystyle\begin{bmatrix}f(x_{1})\\\ \vdots\\\
f(x_{T})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(x_{1})\\\ \vdots\\\
m(x_{T})\end{bmatrix},\begin{bmatrix}k(x_{1},x_{1})\dots k(x_{1},x_{T})\\\
\vdots\\\ k(x_{T},x_{1})\dots k(x_{T},x_{T})\end{bmatrix}\ \right).$ (2)

Here, $\mathcal{N}$ denotes the normal distribution and $T$ is the number of
samples. The input design parameters or samples are denoted as
$x_{1},\ldots,x_{T}$. The function $m:\mathbb{R}^{N}\to\mathbb{R}$ is the
user-defined mean function, often chosen to be $0$ constant function. The
covariance functions are $k:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$.
The posterior Gaussian probability distribution of a new sample point $x$ can
be inferred using Bayes’ rule:

$\displaystyle f(x)|x,x_{1:T},f(x_{1:T})\sim\mathcal{N}(\mu(x),\sigma^{2}(x))$
(3) $\displaystyle\mu(x)\ =\
k(x,x_{1:T})k(x_{1:T},x_{1:T})^{-1}\left(f(x_{1:T})-m(x_{1:T})\right)+m(x)$
$\displaystyle\sigma^{2}(x)\ =\
k(x,x)-k(x,x_{1:T})k(x_{1:T},x_{1:T})^{-1}\sigma(x_{1:T},x)\ ,$

where $x_{1:T}=[x_{1},\dots,x_{T}]$ and
$k(x_{1:T},x_{1:T})=[k(x_{1}),\dots,k(x_{T});\dots;k(x_{T},x_{1}),\dots,k(x_{T},x_{T})]$.
The function $\mu$ and $\sigma^{2}$ are the posterior mean and variance,
respectively. The covariance functions, or the kernels, of the Gaussian
process impact the accuracy of the surrogate model significantly. Choices of
the functions include the squared exponential covariance function, matern
functions, etc. The hyperparameters of the kernel, often denoted as $\theta$,
are obtained through the fitting or training process of the Gaussian process.

Selection of a suitable acquisition function is critical to an efficient and
effective Bayesian optimization algorithm as it determines the location of the
new samples. One of the most popular acquisition functions is expected
improvement (EI) Brochu et al. (2010), which is computed from the conditional
expectation $\mathbb{E}$ of an improvement function defined as

$\displaystyle I(x)=\max\\{f(x^{+})-f(x),0\\},$ (4)

where

$\displaystyle x^{+}=\underset{\begin{subarray}{c}x_{i}\in
x_{1:T}\end{subarray}}{\text{argmin}}f(x_{i}).$ (5)

In other words, $x^{+}$ has the best objective value among existing samples.
The expected improvement is the expectation of (4) conditioned on existing
samples. An important advantage of the expected improvement function is that
is has a closed form:

$\displaystyle
EI(x)=\begin{cases}(f(x^{+})-\mu(x))\Phi(Z)+\sigma(x)\phi(Z),&\sigma(x)>0,\\\
0,&\sigma(x)=0,\end{cases}$ (6)

where

$\displaystyle Z\ =\ \begin{cases}\frac{f(x^{+})-\mu(x)}{\sigma(x)},&\
\sigma(x)>0\\\ 0,&\ \sigma(x)=0\ .\end{cases}$ (7)

A trade-off parameter $\xi\geq 0$ can be added to (6) for more control over
global search and local optimization Lizotte (2008); Jones (2001). The
functions $\phi$ and $\Phi$ are the PDF and CDF of the standard normal
distribution, respectively. Therefore, they are readily available in most of
the coding languages that one might use to implement Bayesian optimization.

The next sample point $x_{k}$ in the Bayesian optimization algorithm is chosen
by maximizing the acquisition function over $C$. At the $k$th iteration, this
means

$\displaystyle x_{k}=\underset{\begin{subarray}{c}x_{i}\in
C\end{subarray}}{\text{argmin}}EI_{k}(x),$ (8)

where $EI_{k}$ is the expected improvement function using samples accumulated
up to the $k$th iteration. In order to solve (8), optimization algorithms
including L-BFGS or random search can be used.

Expected constrained improvement (ECI) Bernardo et al. (2011); Gelbart et al.
(2014); Gardner et al. (2014) is one of the most widely adopted acqusition
functions for constrained optimization problems due to its effectiveness and
ease to implement. The constrained improvement function is defined as

$\displaystyle I_{C}(x)=\Delta(x)\max\\{f(x^{+})-f(x),0\\}=\Delta(x)I(x),$ (9)

where $\Delta(x)\in\\{0,1\\}$ is a feasibility indicator function that has the
value $1$ if $x$ is a feasible point and $0$ otherwise. The optimal existing
sample $x^{+}$ is now

$\displaystyle x^{+}=\underset{\begin{subarray}{c}x_{i}\in x_{1:T}\\\ c(x)\geq
0\end{subarray}}{\text{argmin}}f(x_{i}).$ (10)

It is clear from (10) that a feasible point in the initial data set is
required. The expected value of the constrained improvement function has a
simple closed form

$\displaystyle EI_{C}(x)=PF(x)EI(x),$ (11)

where $PF(\cdot)$ is the probability of a feasible $x$ and it follows a simple
univariate Gaussian CDF. For multiple inequality constraints, i.e., $m>1$,
they are often assumed to be conditionally independent and the probability of
$x$ being feasible is computed as

$\displaystyle p(c(x)\geq 0)=\prod p(c_{i}(x)\geq 0),$ (12)

where $c_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is the $i$th component of the
constraint functions.

##  3 Merit constrained Bayesian optimization

We introduce our proposed acqusition functions and algorithms in this section.
Our idea is to construct the acquisition function based on merit functions
that have been widely adopted in modern nonlinear optimization methods such as
sequential quadratic programming (SQP) Nocedal and Wright (2006). While in
theory all merit functions can be applied here, we focus on the penalty merit
function in this paper.

###  3.1 Expected merit improvement

The penalty merit function can be defined as

$\displaystyle\varphi(x)=f(x)+\alpha c^{+}(x):=f(x)+\alpha\max\\{-c(x),0\\},$
(13)

where $\alpha\geq 0$ is a penalty parameter. It is generally well-understood
how $\alpha$ should be updated Nocedal and Wright (2006). Applying (13) to
define the improvement function has been considered previously Zhao and Xu
(2024); Jiao et al. (2019). However, as is evident in its form, this would
lead to a maximum of non-Gaussian distribution in the form of

$\displaystyle I_{m,0}(x)=\max\\{f(x^{+})+\alpha c^{+}(x^{+})-[f(x)+\alpha
c^{+}(x)],0\\},$ (14)

where

$\displaystyle x^{+}=\underset{\begin{subarray}{c}x_{i}\in
x_{1:T}\end{subarray}}{\text{argmin}}f(x_{i})+\alpha
c^{+}(x_{i})=\underset{\begin{subarray}{c}x_{i}\in
x_{1:T}\end{subarray}}{\text{argmin}}\varphi(x_{i}).$ (15)

The expected value of (14) no longer has a closed form and requires expensive
sampling to estimate. To resolve this issue, we look back at Jones et al.
(1998), where the authors argued for the expected value of an improvement
function, i.e., EI, instead of the expected value of the objective itself for
the much superior globalization property of the former. However, for the
constraint violation functions $c^{+}$, a maximum function similar to that in
the improvement function is already present. Hence, the global search for $c$
could also exist in the expectation of $c^{+}$. Therefore, we propose to relax
form (14) into two merit improvement functions. The first form is

$\displaystyle I_{m,1}(x)=$ $\displaystyle\max\\{f(x^{+})-f(x),0\\}+\alpha
c^{+}(x^{+})-\alpha c^{+}(x)$ (16) $\displaystyle=$
$\displaystyle\max\\{f(x^{+})-f(x),0\\}+\alpha\max\\{-c(x^{+}),0\\}-\alpha\max\\{-c(x),0\\}.$

And the second form is

$\displaystyle I_{m,2}(x)=$ $\displaystyle f(x^{+})-f(x)+\alpha
c^{+}(x)-\alpha c^{+}(x)$ (17) $\displaystyle=$
$\displaystyle\varphi(x^{+})-f(x)-\alpha\max\\{-c(x),0\\}.$

The essential idea of (16) and (17) is that since the expected values of the
constraint violation functions already use the uncertainty information from
$c$ via the maximum functions, we might be able to maintain good globalization
property of the acqusition function while removing the additional maximum
function in (14). We point out that the second form (17) does not use
uncertainty information from the objective, thus could be less effective if
the objective functions have many local minimum.

We call the expected value of (16) the form 1 expected merit improvement
(EMI), denoted as $EI_{M,1}$, which has the closed form

$\displaystyle EI_{M,1}(x)=$
$\displaystyle(f(x^{+})-\mu_{f}(x))\Phi(Z_{f})+\sigma_{f}\phi(Z_{f})+\alpha
c^{+}(x^{+})-\alpha[-\mu_{c}(x)\Phi(Z_{c})+\sigma_{c}\phi(Z_{c})],$ (18)
$\displaystyle=$
$\displaystyle(f(x^{+})-\mu_{f}(x))\Phi(Z_{f})+\sigma_{f}\phi(Z_{f})+\alpha
c^{+}(x^{+})+\alpha\mu_{c}(x)\Phi(Z_{c})-\alpha\sigma_{c}\phi(Z_{c}),$

where

$\displaystyle Z_{f}=\frac{f(x^{+})-\mu_{f}(x)}{\sigma_{f}},\
Z_{c}=\frac{-\mu_{c}(x)}{\sigma_{c}}.$ (19)

The expression for form 2 EMI of (17) is

$\displaystyle EI_{M,2}(x)=$
$\displaystyle\varphi(x^{+})-\mu_{f}(x)-\alpha[-\mu_{c}(x)\Phi(Z_{c})+\sigma_{c}\phi(Z_{c})]$
(20) $\displaystyle=$
$\displaystyle\varphi(x^{+})-\mu_{f}(x)+\alpha\mu_{c}(x)\Phi(Z_{c})-\alpha\sigma_{c}\phi(Z_{c}),$

where $Z_{c}$ is the same as in (19). Both acqusition functions can be easily
extended to multiple constraints. The form 1 EMI becomes

$\displaystyle EI_{M,1}(x)$
$\displaystyle=(f(x^{+})-\mu_{f}(x))\Phi(Z_{f})+\sigma_{f}\phi(Z_{f})+\alpha\sum_{j=1}^{m}c^{+}_{j}(x^{+})$
(21)
$\displaystyle+\alpha\sum_{j=1}^{m}[\mu_{c_{j}}(x)\Phi(Z_{c_{j}})-\sigma_{c_{j}}\phi(Z_{c_{j}})],$

where for $j=1,\dots,m$,

$\displaystyle Z_{c_{j}}=\frac{-\mu_{c_{j}}(x)}{\sigma_{c_{j}}}.$ (22)

Similarly, the second form is

$\displaystyle EI_{M,2}(x)=$
$\displaystyle\varphi(x^{+})-\mu_{f}(x)+\alpha\sum_{j=1}^{m}[\mu_{c_{j}}(x)\Phi(Z_{c_{j}})-\sigma_{c_{j}}\phi(Z_{c_{j}})].$
(23)

The penalty merit function (13) with multiple constraints can be defined as

$\displaystyle\varphi(x)=f(x)+\alpha\sum_{j=1}^{m}\max\\{-c_{j}(x),0\\}.$ (24)

The optimal existing sample is computed using the same formula (15) with the
updated merit function (24). The MCBO algorithms are given in Algorithm 1.

Algorithm 1 Merit constrained Bayesian optimization

1:Choose parameter $\alpha\geq 0$ and its update rule. Choose either (21) or
(23) as the acquisition function. Choose $N_{0}$ initial samples
$x_{i},~{}i\in[0,N_{0}]$. Compute $f(x_{i})$ and $c(x_{i})$.

2:Train Gaussian Process surrogate models for the objective and constraint
functions on the initial samples.

3:for $k=0,1,2,\dots$ do

4: Evaluate the acquisition function and find the new sample point based on
(8).

5: Run (numerical) experiments with the new sample $x_{k}$ if necessary.

6: Evaluate the objective $f(x_{k})$ and constraint functions $c(x_{k})$.
Update $\alpha$ if needed.

7: Retrain the surrogate models with the addition of the new sample $x_{k}$,
$f(x_{k})$ and $c(x_{k})$.

8: Solve the optimization problem (1) with the updated surrogate model.

9: if Stopping criteria satisfied then Exit

We include line $5$ in the algorithm because new numerical simulations might
need to be run with $x_{k}$, as in many engineering applications. Line $8$ of
Algorithm 1 can be solved through exhaustive search on a random, sufficiently
large set of sample points in $C$ or gradient-based methods such as L-BFGS.
The specific stopping criteria and the update rule for $\alpha$ are left to be
chosen by users to allow maximum flexibility of the algorithms. We remark that
the algorithmic parameter $\alpha$ can be changed, e.g., based on whether the
constraints are satisfied, throughout the optimization process. Additionally,
each constraint can have a different $\alpha$ for better performance.

###  3.2 A unified algorithm with with ECI

For the penalty merit function to produce equivalent solution to the original
problem (1), the parameter $\alpha$ needs to be large enough (Nocedal and
Wright, 2006, Theorem 17.3). As a result, the additive structure of the merit
function in (21) and (23) might lead to more exploration within the region
that current constraint surrogate models expect to be feasible. If the
boundary of the feasible regions are complex, the algorithms could become less
efficient. Meanwhile, the multiplication structure in the ECI method, without
the penalty $\alpha$, can lead to more exploration where the current
constraint surrogate models deem close to the boundary but infeasible.
Therefore, we propose a unified acqusition function that combines ECI with EMI
through algorithmic parameters, referred to as the unified expected
constrained improvement (UECI). Its improvement function is given by

$\displaystyle I_{u}(x)=$
$\displaystyle(1-\beta)\Delta(x)\max\\{f(x^{+})-f(x),0\\}$ (25)
$\displaystyle+\beta\left[\max\\{f(x^{+})-f(x),0\\}+\alpha\sum_{j=1}^{m}c_{j}^{+}(x^{+})-\alpha\sum_{j=1}^{m}c^{+}(x)\right],$

where $\beta\in[0,1]$ is a parameter. When $\beta=0$, we recover ECI (9), and
when $\beta=0$, (25) reduces to (16). We note that many other combinations are
possible using, e.g., (22) and other algorithmic parameter design. The UECI
expression is

$\displaystyle EI_{UC}=$ $\displaystyle(1-\beta)PF(x)EI(x)+\beta EI_{m,1}(x),$
(26)

based on (11), (12), and (21). The UECI function continues to be closed-
formed, making its implementation straightforward. In this paper, we devise a
simple update rule for $\beta$ that relies on the number of feasible points.
Specifically, let $n_{f}$ be the number of existing feasible samples and
$N_{f}$ the threshold. A unified constrained Bayesian optimization (UCBO)
algorithm based on UECI (26) is given in Algorithm 2

Algorithm 2 Unified constrained Bayesian optimization

1:Choose parameter $\alpha\geq 0$ and $N_{f}\geq 1$. Set $\beta=1$. Choose
$N_{0}$ initial samples $x_{i},~{}i\in[0,N_{0}]$. Compute $f(x_{i})$ and
$c(x_{i})$.

2:Train Gaussian Process surrogate models for the objective and constraint
functions on the initial samples.

3:for $k=0,1,2,\dots$ do

4: if $n_{f}\geq N_{f}$ then

5: Set $\beta=0$.

6: Evaluate the acquisition function (26) and find the new sample point with
(8).

7: Run (numerical) experiments with the new sample $x_{k}$ if necessary.

8: Evaluate the objective $f(x_{k})$ and constraint functions $c(x_{k})$.
Update $\alpha$ if needed.

9: Retrain the surrogate models with the addition of the new sample $x_{k}$,
$f(x_{k})$ and $c(x_{k})$.

10: Solve the optimization problem (1) with the updated surrogate model.

11: if Stopping criteria satisfied then

12: Exit

Algorithm 2 offers a simplistic switch between ECI and MCBO methods. The
choices for $\beta$ do not need to be binary. We leave more complex designs
for update rules of both $\alpha$ and $\beta$, as well as algorithms, for
future work. Algorithm 1 and 2 have been implemented in Python using the
scikit-learn machine learning library Pedregosa et al. (2011); Buitinck et al.
(2013) and MATLAB.

##  4 Experiments

In this section, we demonstrate and evaluate the performance of MCBO and UCBO
algorithms on three synthetic problems that have been studied and compared
before Gardner et al. (2014); Gramacy et al. (2016); Picheny et al. (2016);
.Ariafar et al. (2019). Additionally, we apply MCBO and UCBO algorithms to a
data-driven ICF design problem based on complex finite element simulation data
and practical design objective and constraints.

For the synthetic problems, we compare our proposed algorithms with the state-
of-the-art constrained Bayesian optimization methods ADMMBO method and the ECI
Gardner et al. (2014) method with additional random sampling. As reported in
.Ariafar et al. (2019), ADMMBO is considered to be one of the best performing
CBO algorithms, as it is built upon previous AL methods. Meanwhile, ECI is one
of the most widely adopted CBO methods thanks to its simplicity and
effectiveness. For a fair comparison, the ADMMBO results in this paper is
obtained using the open-source code the authors of .Ariafar et al. (2019)
published, where the ADMMBO parameters and initial data sets are unchanged to
maximize its performance. Therefore, readers can find the detailed parameter
choices and their discussion for ADMMBO in .Ariafar et al. (2019). If no
initial samples are feasible, we allow the ECI implementation to call Latin
hypercube sampling to generate new samples until a feasible one is found. The
MCBO, UCBO, and ECI algorithms are implemented in both Python and MATLAB,
while ADMMBO is currently available in MATLAB.

Each method is run $100$ times to obtain a more comprehensive comparison of
performances. ADMMBO has two initial samples for the objective and
constraints, respectively, while four initial Latin hypercube samples are used
for other algorithms. The median of the best feasible objectives among the
$100$ runs is reported against the number of iterations. Additionally, we
report the $25$ and $75$ percentiles of the best feasible objectives from our
proposed algorithms to illustrate their robustness. For the ADMMBO method, the
iteration number refers to its inner iteration, which is given a budget of
$100$ by the authors of .Ariafar et al. (2019). Thanks to a stopping
criterion, ADMMBO sometimes converges before the budget is used up. The
budgets of other algorithms are chosen based on when the best feasible
objectives stabilize, as more iterations do not change their performances.

We choose to compare the number of iterations since the evaluation of
objective and constraint functions can be done independently. More
importantly, in many engineering applications, e.g., ICF design optimization,
a new sample in the design space requires a new complex and computationally
expensive numerical simulation, while the design objective and constraints
evaluations are analytical and thus much less costly. We do point out that
ADMMBO has the advantage that only one function, either the objective or the
constraint, is evaluated at one inner iteration.

###  4.1 Example with small feasible region

The mathematical expression for the first example is

$\displaystyle\underset{\begin{subarray}{c}x\end{subarray}\in
C}{\text{minimize}}$ $\displaystyle\sin(x_{1})+x_{2},$ (27) subject to
$\displaystyle-\sin(x_{1})\sin(x_{2})-0.95\geq 0,$

where $C=[0,6]\times[0,6]$. This is a classic problem that has been used in
literature due to its small feasible region and the existence of distinct
local and global minimum. To demonstrate this point, the feasible region and
the objective contour are plotted in figure 1. The local minimum has a value
of $5.4$ and the global minimum has the value of $0.25$, obtained through
exhaustive random sampling of $10000$ points.

Figure 1: Contour plots for the objective function (left) and constraint
function (right) for example 1. The feasible region is marked on the
constraint plot, along with the global optimum ($*$ sign).

Given the small feasible region, we choose $\alpha=20$ for form 1 (21) of
MCBO, $\alpha=5$ for form 2 (23), and $\alpha=20,N_{f}=2$ for UCBO. MCBO and
UCBO algorithms are given $60$ iterations, while ECI and ADMMBO are given
$100$.

Figure 2 shows the median of best feasible objective against number of
iterations among $100$ runs. Form 2 of MCBO and UCBO converge the fastest to
the optimal value, followed closely by form 1 of MCBO, while ECI is the
slowest. This is likely due to the small feasible region that ECI is incapable
of finding effectively. All algorithms succeed in escaping the local minimum
and searching in the neighborhood of the global optimum. The $25$ and $75$
percentiles of the best feasible objectives from MCBO and UCBO algorithms over
the $100$ runs illustrate the robustness and reliability of the proposed
algorithms.

Figure 2: Optimization result of MCBO, UCBO, ECI (median) and ADMMBO (median)
algorithms for example 1, shown in figure 1. The best feasible objective is
plotted against the number of iterations, or equivalently, number of samples.

It is worth pointing out that a more sophisticated update rule for $\alpha$
and $\beta$ is likely to improve performance further for the proposed
algorithms, a topic for future work.

###  4.2 Example with multiple constraints

Next, we consider the test problem

$\displaystyle\underset{\begin{subarray}{c}x\end{subarray}\in
C}{\text{minimize}}$ $\displaystyle x_{1}+x_{2},$ (28) subject to
$\displaystyle 0.5\sin(2\pi(x_{1}^{2}-2x_{2}))+x_{1}+2x_{2}-1.5\geq 0,$
$\displaystyle-x_{1}^{2}-x^{2}_{2}+1.5\geq 0,$

where $C=[0,1]\times[0,1]$. The feasible region and the objective contour are
plotted in figure 3. It is evident from the plot that the direction of a
smaller objective is in agreement with the feasible region of the second
constraint. Additionally, the feasible region is not as restrictive as the one
in example 1.

Figure 3: Contour plots for the objective function (left) and the two
constraint functions (middle and right) for example 2. The feasible region is
marked with black line on the objective contour, along with the global optimum
($*$ sign). The contour curve $c(x)=0$ is in black on each constraint contour
plot as well.

We choose $\alpha=[2,0.02]$ for form 1 of MCBO, $\alpha=[25,25]$ for form 2 of
MCBO, and $\alpha=[100,0.1],N_{f}=1$ for UCBO. ADMMBO is allowed the same
$100$ inner iterations without modifying the original code to ensure fair
comparison. Based on their performances, form 1 of MCBO, UCBO, and ECI are run
for $60$ iterations while form 2 has a budget of $80$ iterations that matches
the iterations reported by ADMMBO upon exit.

The median of the best feasible objectives are plotted in figure 4, together
with the best feasible objectives at $25$ and $75$ percentiles for MCBO and
UCBO algorithms. Form 2 MCBO finds feasible samples and converges to optimal
solution faster than ADMMBO, while form 2 initial outperforms ADMMBO but
eventually find the optimum in similar number of iterations. The $25\%$ and
$75\%$ curves demonstrate the robustness and stability of the two proposed
algorithms. UCBO is the best performing algorithm among all those tested.

Figure 4: Optimization result of MCBO, UCBO, ECI (median), and ADMMBO (median)
algorithms for example 2, shown in figure 3. The best feasible objective is
plotted against the number of iterations.

Perhaps the surprising result is the excellent performance of ECI for this
problem in our experiments, contrary to findings in previous papers .Ariafar
et al. (2019); Gramacy et al. (2016); Hernández-Lobato et al. (2015). For
example, in .Ariafar et al. (2019), ECI does not have feasible samples for the
$100$ runs until $40$ iterations, a significantly worse performance than ours.
In our experiments, given that the feasible region is relative large in this
low-dimensional problem, ECI, with the help of continued Latin hypercube
sampling, is able to have a feasible sample within the first eight samples. It
is then capable of finding global minimum very efficiently. Since it is
unclear precisely how other authors overcome the lack of initial feasible
samples for ECI in their experiments, we contemplate that our approach of
sampling when there is no feasible sample could be the reason for the improved
performance of ECI. The UCBO method can fully take advantage the capability of
ECI to explore outside of feasible region while working with infeasible
initial samples and locating feasible ones quickly.

###  4.3 Four-dimensional problem

To examine the performance of our proposed algorithms for problems in higher
dimensions, we test a four-dimensional problem that previously appeared in
Picheny et al. (2016); .Ariafar et al. (2019)

$\displaystyle\underset{\begin{subarray}{c}x\end{subarray}\in
C}{\text{minimize}}$ $\displaystyle\sum_{i=1}^{4}x_{i},$ (29) subject to
$\displaystyle\frac{1}{0.8387}\left[\sum_{i=1}^{4}E_{i}\exp(-\sum_{j=1}^{4}A_{ji}(x_{j}-P_{ji})^{2})-1.1\right]\geq
0,$

where $C=[0,1]^{4}$ and

$\displaystyle E=\begin{bmatrix}&1\\\ &1.2\\\ &3\\\ &3.2\end{bmatrix},\
A=\begin{bmatrix}&10\ &0.05\ &3\ &17\\\ &3\ &10\ &3.5\ &8\\\ &17\ &17\ &1.7\
&0.05\\\ &3.5\ &0.1\ &10\ &10\end{bmatrix},\ P=\begin{bmatrix}&0.131\ &0.232\
&0.234\ &0.404\\\ &0.169\ &0.413\ &0.145\ &0.882\\\ &0.556\ &0.83\ &0.352\
&0.873\\\ &0.012\ &0.373\ &0.288\ &0.574\\\ \end{bmatrix}.$ (30)

We choose an iteration budget of $50$ for all five algorithms. We let
$\alpha=0$ for the first $10$ iterations of form 1 MCBO and $\alpha=0.01$ for
the remaining 40 iterations, while we set $\alpha=0.01$ for form 2 MCBO.
Similarly, we let $\alpha=0.01$ and $N_{f}=2$ for UCBO. The result of the
optimization is shown in figure 5, where the same curves as in example 1 and 2
are plotted. Among the median results, form 1 MCBO algorithm is the first to
find feasible samples. Moreover, all algorithms quickly find values close to
the global minimum, except for ADMMBO, which might be caused by the constraint
tolerance $0.01$ set by the ADMMBO run. Therefore, the proposed algorithms are
capable of top-tier performance in higher dimensions as well.

Figure 5: Optimization result of MCBO, UCBO, ECI (median) and ADMMBO (median)
algorithms for example 3. The best feasible objective is plotted against the
number of iterations.

It is worth noting that the ADMMBO results for all three synthetic problems
here is very similar to those presented in .Ariafar et al. (2019), where it is
the best performing algorithm among all those tested.

###  4.4 ICF design problem

Finally, we present a data-driven ICF design optimization problem, where the
objective and constraint functions are results from complex simulation of
multiphysics problems with no analytical expressions. The objective and
constraint functions are based on simulation data collected using HYDRA, a
multi-physics simulation code developed at LLNL Marinak et al. (2001). To
perform an efficient study, HYDRA simulation is not called directly by our
algorithms. Rather, a pre-trained neural networks fitted to a HYDRA simulation
database of one-dimensional capsule simulations of design variations around
NIF shot N210808 is used. It is the first ICF experiment to exceed the Lawson
fusion condition and generate more than 1 MJ of energy Abu-Shawareb et al.
(2022). The optimization process would be prohibitively expensive if HYDRA is
run natively.

We select design variables ‘t_2nd’ and ‘sc_peak’, which are the timing of the
second shock and the strength of the peak of the radiation drive. For a figure
that illustrates the relationship between the two variables and the radiation
temperature, see Wang et al. (2024). The design objective is to maximize the
nuclear yield, or equivalently minimizing the negative yield. The two
constraint functions are chosen to be $c_{1}$ = ‘adiabat’ and $c_{2}$= ‘vImp’.
The implosion velocity vImp (km/sec) is the maximum inward velocity of the
fusion fuel achieved during the implosion, and the adiabat is a unit-less
measure $>1$ of the entropy of the fuel at the time of peak implosion
velocity. In one-dimensional models (like simulated here), implosions with
lower values of adiabat and higher values of velocity tend to achieve higher
yield; however, more complex models show that additional performance-limiting
physics can turn on if either the adiabat becomes too low or the velocity
becomes too high. As such, it is often desirable to search for designs that
have constraints on these functions. Therefore, we impose the constraints
$c_{1}(x)\geq 4.25$ and $c_{2}(x)\leq 350$.

The contour plots for the objective (negative yield) and constraints are given
in figure 6, along with the feasible region. The optimal feasible objective
value is obtained through exhaustive sampling of the two-dimensional design
space. As evident in figure 6, the feasible region is relatively small for
this problem.

Figure 6: Contour plots for the objective function (left) and the two
constraint functions (middle and right) from HYDRA data set. The feasible
region is marked with black dashed lines on the objective contour, along with
the global optimum ($*$ sign). The contour curve $c(x)=0$ is in black curve on
each constraint contour plot as well.

We adopt to only apply MCBO, UCBO, and ECI algorithms for this problem since
they are the most straightforward to implement in an existing and complex code
base. For all three proposed algorithms, we choose $\alpha=[1,1]$ and set
$N_{f}=3$ for UCBO due to the small feasible region. We allow ECI to randomly
sample until a feasible point is found. A total budget of $50$ iterations is
assigned to all algorithms tested. The result is given in figure 7, where the
MCBO and UCBO algorithms behave similarly but clearly outperform ECI thanks to
its ability to make progress on constraint violation from the beginning.

Figure 7: Optimization result of two MCBO, UCBO and ECI (median) algorithms
for the ICF design problem in figure 6. The mean/25%/75% best feasible
objectives of MCBO and UCBO methods are plotted against the number of
iterations.

The model ICF design problem presented here can be extended to a full problem
with eight or as high as thirty dimensions. An eight-dimensional constrained
ICF design problems is currently being developed.

##  5 Conclusion

In this paper, we propose constrained Bayesian optimization methods with merit
expected improvement acquisition functions that can resolve a large number of
issues that constrained expected improvement encounter. Through relaxation on
the improvement function of the constraint violation functions, we can
maintain closed-form expressions for the new acquisition functions, while
making progress at infeasible samples. The additive structure between the
objective and constraint violation functions allows expectation to be taken
efficiently on conditionally dependent constraint functions. Moreover, the
proposed algorithms do not require any new routines to implement compared to
unconstrained Bayesian optimization.

We further propose a unified expected constrained improvement function and
design a constrained Bayesian optimization algorithm that can take advantage
of the benefits of expected constraint improvement function if needed. The
unified approach can be seen as a framework where other methods can be
integrated.

Some questions that have not been answered in this paper include the
algorithmic parameter update rule for $\alpha$ and $\beta$. More numerical
experiments in higher dimension would further validate the performance of the
proposed algorithms. The authors will work on these topics in the future.

## Acknowledgments and Disclosure of Funding

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under contract DE–AC52–07NA27344 and
the LLNL-LDRD Program under Project tracking No. 23-ERD-017. Release number
LLNL-CONF-861915-DRAFT. The data and codes that support the findings of this
study are available from the corresponding author upon request.

## References

  * Wang et al. [2023] X. Wang, Y. Jin, S. Schmitt, and M. Olhofer.  Recent advances in bayesian optimization.  _ACM Comput. Surv._ , 55(13s), jul 2023.
  * Wang and Zabaras [2004] J. Wang and N. Zabaras.  A Bayesian inference approach to the inverse heat conduction problem.  _International Journal of Heat and Mass Transfer_ , 47:3927–3941, August 2004.
  * Mathern et al. [2021] A. Mathern, O. S. Steinholtz, A. Sjöberg, et al.  Multi-objective constrained Bayesian optimization for structural design.  _Structural and Multidisciplinary Optimization_ , 63:689–701, February 2021.
  * Wang and Papadopoulos [2023] J. Wang and P. Papadopoulos.  Optimization of process parameters in additive manufacturing based on the finite element method.  _arXiv preprint arXiv:2310.15525_ , 2023.
  * Calandra et al. [2016] R. Calandra, A. Seyfarth, J. Peters, and M. P. Deisenroth.  Bayesian optimization for learning gaits under uncertainty.  _Annals of Mathematics and Artificial Intelligence_ , 76:5–23, February 2016.
  * Wang et al. [2024] J. Wang, N. Y. Chiang, A. Gillette, and J. L. Peterson.  A multifidelity Bayesian optimization method for inertial confinement fusion design.  _Phys Plasmas (in press)_ , 2024.
  * Zuluaga et al. [2013] M. Zuluaga, G. Sergent, A. Krause, and M. Püschel.  Active learning for multi-objective optimization.  In _International Conference on Machine Learning_ , pages 462–470, February 2013.
  * Bernardo et al. [2011] J. Bernardo, M. J. Bayarri, J. O. Berger, A. P. Dawid, D. Heckerman, A. F. Smith, and M. West.  Optimization under unknown constraints.  _Bayesian Statistics_ , 9:229, October 2011.
  * Schonlau et al. [1998] M. Schonlau, W. J. Welch, and D. R. Jones.  Global versus local search in constrained optimization of computer models.  _Lecture notes-monograph series_ , pages 11–25, 1998.
  * Gelbart et al. [2014] M. A. Gelbart, J. Snoek, and R. P. Adams.  Bayesian optimization with unknown constraints.  _arXiv preprint arXiv:1403.5607_ , 2014.
  * Gardner et al. [2014] J. R. Gardner, M. J. Kusner, Z. Xu, K. Q. Weinberger, and J. P. Cunningham.  Bayesian optimization with inequality constraints.  In _Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32_ , ICML’14, pages II–937–II–945. JMLR.org, 2014.
  * Letham et al. [2019] B. Letham, B. Karrer, G. Ottoni, and E. Bakshy.  Constrained Bayesian optimization with noisy experiments.  2019\.
  * Eriksson and Poloczek [2021] D. Eriksson and M. Poloczek.  Scalable constrained Bayesian optimization.  In _International Conference on Artificial Intelligence and Statistics_ , pages 730–738. PMLR, 2021.
  * Picheny [2014] V. Picheny.  A stepwise uncertainty reduction approach to constrained global optimization.  In _Artificial intelligence and statistics_ , pages 787–795. PMLR, 2014.
  * Hernández-Lobato et al. [2015] J. M. Hernández-Lobato, M. Gelbart, M. Hoffman, R. Adams, and Z. Ghahramani.  Predictive entropy search for Bayesian optimization with unknown constraints.  In _International conference on machine learning_ , pages 1699–1707. PMLR, 2015.
  * Zhang et al. [2023] F. Zhang, Z. Zhu, and Y. Chen.  Constrained Bayesian optimization with adaptive active learning of unknown constraints.  _arXiv preprint arXiv:2310.08751_ , 2023.
  * Gramacy et al. [2016] R. B. Gramacy, G. A. Gray, S. L. Digabel, H. K. Lee, P. R. Ranjan, G. Wells, and S. M. Wild.  Modeling an augmented lagrangian for blackbox constrained optimization.  _Technometrics_ , 58(1):1–11, 2016.
  * Picheny et al. [2016] V. Picheny, R. B. Gramacy, S. Wild, and S. L. Digabel.  Bayesian optimization under mixed constraints with a slack-variable augmented lagrangian.  _Advances in neural information processing systems_ , 29, 2016.
  * .Ariafar et al. [2019] S .Ariafar, J. Coll-Font, D. Brooks, and J. Dy.  ADMMBO: Bayesian optimization with unknown constraints using ADMM.  _Journal of Machine Learning Research_ , 20(123):1–26, 2019.
  * Jiao et al. [2019] R. Jiao, S. Zeng, C. Li, Y. Jiang, and Y. Jin.  A complete expected improvement criterion for gaussian process assisted highly constrained expensive optimization.  _Information Sciences_ , 471:80–96, 2019.
  * Zhao and Xu [2024] J. Zhao and J. Xu.  Bayesian optimization via exact penalty.  _Technometrics_ , (just-accepted):1–41, 2024\.
  * Brochu et al. [2010] E. Brochu, V. M. Cora, and N. De Freitas.  A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.  _arXiv preprint arXiv:1012.2599_ , December 2010.
  * Lizotte [2008] D. J. Lizotte.  _Practical bayesian optimization_.  PhD thesis, University of Alberta, Edmonton, Alberta, Canada, 2008.
  * Jones [2001] D. R. Jones.  A taxonomy of global optimization methods based on response surfaces.  _Journal of global optimization_ , 21:345–383, 2001.
  * Nocedal and Wright [2006] J. Nocedal and S. Wright.  _Numerical Optimization_.  Springer-Verlag, New York, 2006.
  * Jones et al. [1998] D. R. Jones, M. Schonlau, and W. J. Welch.  Efficient global optimization of expensive black-box functions.  _Journal of Global optimization_ , 13:455–492, 1998.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Buitinck et al. [2013] L. Buitinck, G. Louppe, M. Blondel, F. Pedregosa, A. Mueller, O. Grisel, V. Niculae, P. Prettenhofer, A. Gramfort, J. Grobler, R. Layton, J. VanderPlas, A. Joly, B. Holt, and G. Varoquaux.  API design for machine learning software: experiences from the scikit-learn project.  In _ECML PKDD Workshop: Languages for Data Mining and Machine Learning_ , pages 108–122, 2013.
  * Marinak et al. [2001] M. M. Marinak, G. D. Kerbel, N. A. Gentile, O. Jones, D. Munro, S. Pollaine, T. R. Dittrich, and S. W. Haan.  Three-dimensional hydra simulations of national ignition facility targets.  _Physics of Plasmas_ , 8(5):2275–2280, 2001.
  * Abu-Shawareb et al. [2022] H. Abu-Shawareb et al.  Lawson criterion for ignition exceeded in an inertial fusion experiment.  _Phys. Rev. Lett._ , 129:075001, Aug 2022.  doi: 10.1103/PhysRevLett.129.075001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.129.075001.