Coupling and selecting constraints in Bayesian optimization under uncertainties

Abstract

We consider Reliability-based Robust Design Optimization (RRDO) where it is sought to optimize the mean of an objective function while satisfying constraints in probability. The high computational cost of the simulations underlying the objective and constraints strongly limits the number of evaluations and makes this type of problems particularly challenging. The numerical cost issue and the parametric uncertainties have been addressed with Bayesian optimization algorithms which leverage Gaussian processes of the objective and constraint functions. Current Bayesian optimization algorithms call the objective and constraint functions simultaneously at each iteration. This is often not necessary and overlooks calculation savings opportunity. This article proposes a new efficient RRDO Bayesian optimization algorithm that optimally selects for evaluation, not only the usual design variables, but also one or several constraints along with the uncertain parameters. The algorithm relies on a multi-output Gaussian model of the constraints. The coupling of constraints and their separated selection are gradually implemented in three algorithm variants which are compared to a reference Bayesian approach. The results are promising in terms of convergence speed, accuracy and stability as observed on a two, a four and a 27-dimensional problem.

Appendices

Appendix A: Main notations and acronyms

Appendix B: Covariance kernel for a nominal input: the hypersphere decomposition.

The hypersphere decomposition (Zhou et al. 2011) is a possible choice to parameterize a discrete covariance kernel. The underlying idea is to map each of the l levels $\{z_1,..,z_p,...,z_l\}$ of the considered discrete variable onto a distinct point on the surface of a l-dimensional hypersphere by:

$$\begin{aligned}{} & {} \phi : \{z_1,..,z_p,...,z_l\}\rightarrow \mathbb {R}^l \\{} & {} \phi (z_p) = \sigma _z [ b_{p,1}, \dots , b_{p,l}]^\top \nonumber \quad \text{ for } p=1,\ldots , l \end{aligned}$$

where $b_{p,r}$ represents the r-th coordinate of the p-th discrete level mapping, and is calculated as follows:

$$\begin{aligned}&b_{1,1} = 1&\\&b_{p,r} = \cos \theta _{p,r} \prod _{k = 1}^{r-1} \sin \theta _{p,k}&\text{ for } r = 1,\dots ,p-1 \\&b_{p,p} =\prod _{k = 1}^{p-1} \sin \theta _{p,k}&\text{ for } p \ne 1 \\&b_{p,r} = 0&\text{ for } r \ge p \ne 1 \end{aligned}$$

with $-\pi \le \theta _{p,r} \le \pi$. It can be noticed that in the equations above, some of the mapping coordinates are arbitrarily set to 0. This allows to avoid rotation indeterminacies (i.e., an infinite number of hyperparameter sets characterizing the same covariance matrix), while also reducing the number of parameters required to define the mapping. The resulting kernel is then computed as the Euclidean scalar product between the hypersphere mappings,

$$\begin{aligned} k(z,z') =\phi (z)^\top \phi (z')~. \end{aligned}$$

The discrete kernel can then be characterized as an $l\times l$ symmetric positive definite matrix $\textbf{T}$ containing the covariance values between the discrete variable levels computed as:

$$\begin{aligned} \textbf{T} = \textbf{L}\textbf{L}^\top \end{aligned}$$

where each element of $\textbf{L}_{i,j}$ is computed as $b_{i,j}$,

$$\begin{aligned} \textbf{L} = \sigma _z \left[ \begin{array}{c c c c c} 1 &{} 0 &{} \dots &{} \dots &{} 0 \\ \cos \theta _{2,1} &{} \sin \theta _{2,1} &{} 0 &{} \dots &{} \dots \\ \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots \\ \cos \theta _{l,1} &{} \sin \theta _{l,1} \cos \theta _{l,2} &{} \dots &{} \cos \theta _{l,l-1} \prod _{r = 1}^{l-2} \sin \theta _{l,r} &{} \prod _{r = 1}^{l-1} \sin \theta _{l,r} \end{array} \right] ~. \end{aligned}$$

In the output-as-input model of the constraints, Equation (3), the matrix $\textbf{T}$ contains the covariance terms related to the constraint index, $k_p(i,j) = \textbf{T}_{ij}~,~i,j = \{1,\ldots ,l\}$.

Appendix C: Probability of feasibility with dependent constraints

The probability of satisfying the coupled constraints intervenes in Eq. (6) for the EFI acquisition criterion. It can be estimated with GPs as:

$$\begin{aligned} \mathbb {P}(C^{(t)}(\textbf{x}) \le 0) \approx \frac{1}{N} \sum \limits _{k=1}^{N} \mathbb {1}_{\big (1 - \alpha - \frac{1}{M} \sum \limits _{j=1}^{M}\mathbb {1}_{\big (\textbf{G}^{(t)}(\textbf{x},\mathbf {u_j},\omega _k) \le \textbf{0} \big )} \le 0\big )}. \end{aligned}$$

(17)

In practice, a set of M instances of the random variables $\textbf{U}$ is sampled from $\rho _U$. Subsequently, N independent multi-output trajectories of $\textbf{G}(\textbf{x},\cdot )$ are simulated at the aforementioned sampled uncertain parameters. The probability of feasibility can finally be computed by simply counting the number of trajectories for which the ratio of samples associated to feasible constraints is larger than $1-\alpha$. Note that a multi-output GP prediction of the constraint vector is defined as feasible when all of its components are below their specific threshold value (which is 0 in this work).

In (7), the improvement is computed with respect to the incumbent best feasible solution $z_{min}^{feas}$. However, the objective function mean is not observed therefore $z_{min}^{feas}$ cannot be read in the data set used to condition the GPs F and $\textbf{G}$. For this reason, in El Amri et al. (2021) the incumbent best feasible solution is defined by taking into account the mean of the GP $Z^{(t)}(\textbf{x})$ and the expected value of the process $C^{(t)}(\textbf{x})$:

$$\begin{aligned} z_{min}^{feas} = \min _{\textbf{x}} m^{(t)}_Z(\textbf{x}) ~\text{ s.t. }~ \mathbb E[ C^{(t)}(\textbf{x}) ] \le 0~. \end{aligned}$$

(18)

Given that the Fubini condition holds since the value of $C(\cdot )$ is bounded by definition, the expected value of the process C can be written

$$\begin{aligned} \mathbb E[C^{(t)}(\textbf{x})]= & {} 1 - \alpha - \mathbb E_\textbf{U}[\mathbb E[1_{\textbf{G} (\textbf{x}, \textbf{U}) \le \textbf{0} } ]] \nonumber \\= & {} 1 - \alpha - \int _{\mathbb {R}^m} \mathbb P(\textbf{G} (\textbf{x}, \textbf{u}) \le \textbf{0}))\rho _\textbf{U} (\textbf{u}) d\textbf{u}\nonumber \\= & {} 1 - \alpha - \int _{\mathbb {R}^m} \Phi \left( \textbf{0} -\textbf{m}_\textbf{G}^{(t)}(\textbf{x},\textbf{u}),\textbf{K}_\textbf{G}^{(t)}(\textbf{x},\textbf{u}) \right) \rho _\textbf{U} (\textbf{u}) d\textbf{u} \end{aligned}$$

where $\Phi (\cdot )$ is the cumulative distribution function of a multivariate Gaussian distribution which, like the univariate version, is estimated numerically (Genz and Bretz 2009). The above Equation and Eq. (18) are used to compute the current feasible minimum. If Eq. (18) yields no solution, $m^{(t)}_Z(\textbf{x})$ at the $\textbf{x}$ providing the largest probability of feasibility is taken as the incumbent optimal solution.

Appendix D: Estimating the proxy of the one-step-ahead feasible improvement variance

During the optimization process, the value of the coordinates $\textbf{u}^{t+1}$ of the point to be added to the training data set $\mathcal {D}$ is computed by minimizing the proxy of the one-step-ahead variance of the EFI at $\textbf{x}_{targ}$ (see Eq. (8)), which is defined as:

$$\begin{aligned}{} & {} S(\textbf{x}_{targ},\tilde{\textbf{u}}) = Var\big ( I^{(t+1)}(\textbf{x}_{targ})\big ) \int _{\mathbb {R}^m} Var\big ( 1_{\{\textbf{G}^{(t+1)}(\textbf{x}_{targ},\textbf{u}) \le \textbf{0} \}} \big ) \rho _{\textbf{U}}(\textbf{u}) d\textbf{u},\\{} & {} \quad =~Var(\big (z_{\min }^{\text{ feas }} - Z^{(t+1)}(\textbf{x}_{targ})\big )^+)\int _{\mathbb {R}^m} Var\big ( 1_{\{\textbf{G}^{(t+1)}(\textbf{x}_{targ},\textbf{u}) \le \textbf{0} \}} \big ) \rho _{\textbf{U}}(\textbf{u}) d\textbf{u}. \nonumber \end{aligned}$$

(19)

In this appendix, we recall some details about the computation of the first term in the previous equation. An expression of the improvement variance $Var\left( I^{(t+1)}(\textbf{x}_{targ})\right)$ which bears some similarities to the one of expected improvement is given in El Amri et al. (2021) and can be expressed in terms of probability and density functions of Gaussian distributions:

$$\begin{aligned} EI^{(s)}(\textbf{x})= & {} \big (z_{\min }^\text {feas} - m_Z^{(s)}(\textbf{x}) \big ) \phi \bigg (\frac{z_{\min }^\text {feas} - m_Z^{(s)}(\textbf{x})}{\sigma _Z^{(s)}(\textbf{x})}\bigg ) + \sigma _Z^{(s)}(\textbf{x}) \Phi \bigg (\frac{z_{\min }^\text {feas} - m_Z^{(s)}(\textbf{x})}{\sigma _Z^{(s)}(\textbf{x})}\bigg ),\\= & {} \psi _{EI}(m_Z^{(s)}(\textbf{x}), \sigma _Z^{(s)}(\textbf{x})). \\ Var\left( I^{(s)}(\textbf{x})\right)= & {} EI^{(s)}(\textbf{x}) \big (z_{\min }^\text {feas} - m_Z^{(s)}(\textbf{x}) - EI^{(s)}(\textbf{x})\big ) \\{} & {} + (\sigma _Z^{(s)})^2(\textbf{x}) \Phi \bigg (\frac{z_{\min }^\text {feas} - m_Z^{(s)}(\textbf{x})}{\sigma _Z^{(s)}(\textbf{x})}\bigg ),\\= & {} \psi _{VI}(m_Z^{(s)}(\textbf{x}), \sigma _Z^{(s)}(\textbf{x})). \end{aligned}$$

At step $t+1$, the training data set $\mathcal {D}^{(t)}$ is enriched by $(\tilde{\textbf{x}}, \tilde{\textbf{u}})$ on which the output $f(\tilde{\textbf{x}},\tilde{\textbf{u}})$ is unknown and represented by $F(\tilde{\textbf{x}},\tilde{\textbf{u}})$.

As a consequence, $Var\left( I^{(t+1)}(\textbf{x}_{targ})\right)$ cannot directly be computed because of the randomness of $m_Z^{(t+1)}(\textbf{x}_{targ})$. Indeed, $m_Z^{(t+1)}(\textbf{x}_{targ})$ is equal to

$E\left( Z(\textbf{x}) \vert F(\mathcal {D}^{(t)})=f^{(t)}, F(\tilde{\textbf{x}}, \tilde{\textbf{u}}) \right)$ and follows:

$$\begin{aligned} m_Z^{(t+1)}(\textbf{x}_{targ}) \sim \mathcal {N}\Bigg ( m_Z^{(t)}(\textbf{x}_{targ}) , \bigg (\frac{\int _{\mathbb {R}^m} k_{*F}^{(t)}(\textbf{x}_{targ},\textbf{u};\tilde{\textbf{x}},\tilde{\textbf{u}}) \rho _\textbf{U}(\textbf{u}) d\textbf{u} }{\sqrt{k_{*F}^{(t)}(\tilde{\textbf{x}},\tilde{\textbf{u}};\tilde{\textbf{x}},\tilde{\textbf{u}})}}\bigg )^2 \Bigg ). \end{aligned}$$

We can note that $\sigma _Z^{(t+1)}(\textbf{x}_{targ})$ is not random as it depends only on the location $(\tilde{\textbf{x}},\tilde{\textbf{u}})$ and not on the function evaluation at this point. By applying the law of total variance, it can be shown that:

$$\begin{aligned} Var \left( I^{t+1}(\textbf{x}_{targ})\right)= & {} \mathbb E\left[ Var\left( \big (z_{\min }^{\text {feas}} - Z(\textbf{x}_{targ})\big )^+ \vert m_Z^{(t+1)}(\textbf{x}_{targ}) \right) \right] \\{} & {} +Var\left[ \mathbb E\left( \big (z_{\min }^{\text {feas}} - Z(\textbf{x}_{targ})\big )^+ \vert m_Z^{(t+1)}(\textbf{x}_{targ}) \right) \right] . \\= & {} \mathbb E\left[ \psi _{VI}(m_Z^{(t+1)}(\textbf{x}), \sigma _Z^{(t+1)}(\textbf{x})) \right] \\{} & {} +Var\left[ \psi _{EI}(m_Z^{(t+1)}(\textbf{x}), \sigma _Z^{(t+1)}(\textbf{x})) \right] . \\ \end{aligned}$$

This calculation is performed numerically using samples of $m_Z^{(t+1)}(\textbf{x}_{targ})$. For the sake of clarity, the reader is referred to the previous work (El Amri et al. 2021) for details regarding the implementation.

Appendix E: Optimization bounds of the industrial test case

The bounds of the optimization problem detailed in Sect. 6.4 are provided in the table below.

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$	$x_8$	$x_9$	$x_{10}$
Lower bound	0.05	0.05	0.05	0.05	0.09	0.09	0.09	0.09	0.45	0.45
Upper bound	0.06	0.06	0.06	0.06	0.11	0.11	0.11	0.11	0.55	0.55

	$x_{11}$	$x_{12}$	$x_{13}$	$x_{14}$	$x_{15}$	$x_{16}$	$x_{17}$	$x_{18}$	$x_{19}$	$x_{20}$
Lower bound	0.45	0.45	30	39	47	56	0.2	0.2	0.2	0.2
Upper bound	0.55	0.55	42	51	59	68	8.0	8.0	8.0	8.0

	$u_1$	$u_2$	$u_3$	$u_4$	$u_5$	$u_6$	$u_7$
Lower bound	0.05	0.0000004	0.98	-5.763	-0.5	19.78	16845
Upper bound	0.80	0.0000012	1.02	5.763	0.5	20.59	17532

REF algorithm

SMCS algorithm

MMCU algorithm

MMCS algorithm

Appendix F: Flow charts of the algorithms

This appendix contains the flowcharts of the four algorithms which are implemented and compared in the body of the article. The REF, SMCS, MMCU and MMCS methods are described in Algorithms 2, 3, 4 and 5, respectively. REF stands for reference and was initially proposed in El Amri et al. (2021). It is a method where the GPs of the constraints and objective function are independent, and the same pair $(\textbf{x}^{t}, \textbf{u}^{t})$ is added to every GP at each iteration. SMCS, which stands for Single Models of the constraints and Constraint Selection, has independent GPs, like the REF algorithm, but only one constraint is updated at each iteration. The random parameters of the objective function and the selected constraint, $\textbf{u}_f$ and $\textbf{u}_g$, are different. MMCU means Multiple Model of the constraints and Common $\textbf{u}$. The MMCU algorithm has a joint model of all the constraints and the same iterate $(\textbf{x}^{t}, \textbf{u}^{t})$ enriches all GPs. Finally, MMCS is the acronym of Multiple Model of the constraints and Constraint Selection. The MMCS algorithm relies on a joint model of the constraints and identifies at each iteration a single constraint and the associated random sample, $\textbf{u}_g$, to carry out the next evaluation.