Efficient Tuning of an Isotope Separation Online System Through Safe Bayesian Optimization with Simulation-Informed Gaussian Process for the Constraints

Abstract

Optimizing process outcomes by tuning parameters through an automated system is common in industry. Ideally, this optimization is performed as efficiently as possible, using the minimum number of steps to achieve an optimal configuration. However, care must often be taken to ensure that, in pursuing the optimal solution, the process does not enter an “unsafe” state (for the process itself or its surroundings). Safe Bayesian optimization is a viable method in such contexts, as it guarantees constraint fulfillment during the optimization process, ensuring the system remains safe. This method assumes the constraints are real-valued and continuous functions. However, in some cases, the constraints are binary (true/false) or classification-based (safe/unsafe), limiting the direct application of safe Bayesian optimization. Therefore, a slight modification of safe Bayesian optimization allows for applying the method using a probabilistic classifier for learning classification constraints. However, violation of constraints may occur during the optimization process, as the theoretical guarantees of safe Bayesian optimization do not apply to discontinuous functions. This paper addresses this limitation by introducing an enhanced version of safe Bayesian optimization incorporating a simulation-informed Gaussian process (GP) for handling classification constraints. The simulation-informed GP transforms the classification constraint into a piece-wise function, enabling the application of safe Bayesian optimization. We applied this approach to optimize the parameters of a computational model for the isotope separator online (ISOL) at the MYRRHA facility (Multipurpose Hybrid Research Reactor for High-Tech Applications). The results revealed a significant reduction in constraint violations—approximately $80\%$—compared to safe Bayesian optimization methods that directly learn the classification constraints using Laplace approximation and expectation propagation. The sensitivity to the accuracy of the simulation model was analyzed to determine the extent to which it is advantageous to use the proposed method. These findings suggest that incorporating available information into the optimization process is valuable for reducing the number of unsafe outcomes in constrained optimization scenarios.

1. Introduction

Expensive objective function evaluations and large input parameter spaces pose challenges in optimizing certain processes [1], such as polymer design [2] and accelerator tuning [3,4,5,6,7,8]. These fields demand efficient optimization algorithms to identify optimal solutions with limited objective function evaluations. Bayesian optimization (BO) is commonly employed for such problems because of its simplicity and effectiveness in identifying optimal parameters with few evaluations [9,10]. BO operates through two main steps. First, it constructs a surrogate probabilistic model of the objective function, typically using a Gaussian process (GP). Second, it suggests the next set of design parameters by maximizing an acquisition function, which balances exploration of the design space and exploitation of current model knowledge. The acquisition function achieves this by utilizing the surrogate model’s predictions and associated uncertainties. Furthermore, in many industrial applications, the optimization process must avoid unsafe conditions, such as excessive contact force in medical robot manipulators [11], collisions that could damage robotic components [12], uncomfortable and possible side effects in deep brain stimulation [13], or harm to particle accelerators during parameter tuning [14].

Avoiding unsafe scenarios in the optimization process requires adherence to safety constraints. Constrained Bayesian optimization addresses this by employing a surrogate probabilistic model to learn these constraints [15,16] and then calculating the probability of satisfying them to select a potentially safe parameter set. However, as the model updates after querying new data, consistently safe outcomes are not guaranteed, especially at the start of the optimization, where only some observations of the constraint are available. Usually, safety constraints are real-valued functions, allowing for quantifying how feasible a query input location is. However, some applications involve classification constraints, where it is only possible to determine whether the constraint belongs to the safe or unsafe class. The challenge of solving optimization problems with classification constraints has been addressed using constrained BO with a probabilistic classifier [17], which calculates the probability of belonging to the safe or unsafe class. The max-value entropy search is an alternative acquisition function that works with classification constraints [18]. Tfaily et al. [17] compared various probabilistic classifiers, including k-nearest neighbor [19], decision trees [20], support vector machines (SVMs) [21,22], and logistic regression [22], within the constrained BO framework. The authors in [23] proposed using probabilistic support vector machines (PSVMs) to model classification constraints, extending SVMs to provide the necessary probability of feasibility. In contrast, Gaussian process classifiers inherently offer classification uncertainty, thus providing the probability of feasibility. Among the variants of GP classifiers, the commonly used alternatives are the Laplace approximation, expectation propagation, and variational GP [24]. Moreover, GP classifiers perform better than other classifiers regarding class prediction error [25,26]. Consequently, in this paper, we use Gaussian process classifiers to model the classification constraints.

Safe Bayesian optimization (SafeOpt) offers a more suitable approach for optimizing applications that must adhere to safety constraints [27] than previously introduced constrained Bayesian optimization methods. SafeOpt uses frequentist uncertainty bounds to form a safe set, from which it selects parameter combinations, progressively expanding it until the entire reachable safe space is covered. This method has been successfully applied in various contexts, including particle accelerators [3,14,28,29], industrial robots [12], and quadrotor control [30]. SafeOpt provides theoretical guarantees of finding an optimal solution within safety constraints, assuming at least one initial feasible solution and that the objective and constraint functions are continuous and bounded within the reproducing kernel Hilbert space (RKHS) [27,31]. Another approach involves using barrier functions to define the safe set, ensuring high probability of constraint fulfillment [32,33,34,35]. A recent SafeOpt enhancement introduced by Rothfuss et al. [36] aims to accelerate algorithm convergence by using meta-learned Gaussian process priors based on existing information. However, these methods depend on strong assumptions about the constraints, which can be difficult to verify in practical applications. Recent advancements have led to alternative SafeOpt formulations that do not rely on RKHS bounds. For instance, Fiedler et al. [37] assume Lipschitz continuity in the objective function and bounded noise to maintain safety guarantees. The approach in [38] adjusts exploration levels based on online conformal predictions of a limited violation rate, thereby avoiding strong assumptions about the constraint functions. These SafeOpt variants have proven effective in optimizing design parameters while ensuring constraint fulfillment, provided the functions are continuous.

However, as classification constraints are not continuous, the theoretical guarantees of SafeOpt do not apply. A relevant example is the isotope separator online (ISOL) technique, which produces pure radioisotopes for nuclear, atomic, astro/biophysics, and medical research. The quality (purity) and quantity (intensity) of the supplied radioisotope—in the form of a radioactive ion beam, commonly abbreviated RIB—depend on accurately tuning the underlying process parameters and their interactions. In operational ISOL facilities like ISOLDE at CERN [39] and the ISAC facility at TRIUMF [40], experienced operators set initial design parameters to produce the requested RIB. A significant challenge is the time-intensive nature of parameter tuning, owing to the numerous design parameters that need adjustment. Additionally, safety constraints must be addressed during tuning to prevent contamination of the ISOL’s transport beamline components with radioactive isotopes. Within the ISOL’s transport beamline, assuming that only certain diagnostic elements are available for measurements, it is only possible to determine whether the beam has arrived at the end location or has been lost within any beamline component, leading to classification constraints.

We propose two approaches to extend SafeOpt’s applicability to classification constraints. The first approach directly learns the true classification constraint using a GP classifier, such as Laplace approximation or expectation propagation [41], as shown in Figure 1a. This method is equivalent to using SafeOpt but for cases where the constraints are not real-valued functions. However, since classification constraints are discontinuous, an adaptation of SafeOpt is introduced to ensure the method remains applicable. The second approach leverages a simulation model of the ground truth system, as illustrated in Figure 1b. The simulation model enables the quantification of constraints and aids in constructing a surrogate model, referred to as a simulation-informed Gaussian process (SI-GP) for the constraints. This method integrates observations from both the true classification constraint and the simulation constraint [42], resulting in a piece-wise function represented by the “SI constraint” block in Figure 1b. By fitting a GP to these piece-wise functions, the simulation-informed GP models facilitate the application of SafeOpt. It is worth noting that both approaches use a conventional GP to model the true objective function, as it is a real-valued function. In this paper, we validate the proposed method for optimizing the future transport beamline of the ISOL@MYRRHA facility at the Belgian Nuclear Research Centre (SCK CEN) [43,44]. As the facility is not yet operational, we rely on a computational model as the ground truth model. Although the SI-GP approach is applied to the transport beamline of ISOL@MYRRHA, it can be extended to any optimization problem involving a simulation model of the constraint. The core concept involves formulating a virtual constraint as a piece-wise function, replacing the discontinuous classification constraint.

As a remark, the methods are not tested on benchmark functions, as they build upon SafeOpt, whose convergence and performance in handling safety constraints have already been demonstrated [14,28]. Instead, this paper aims to show how the method using the SI-GP can enhance safety and reduce constraint violations compared to methods that directly learn classification constraints. This is achieved by comparing safe Bayesian optimization with simulation-informed GP (SI-SafeBO) against two other formulations: one using Laplace approximation and one using expectation propagation to approximate true classification constraints. The evaluation metrics include convergence speed and the number of constraint violations. The contributions of this work are summarized as follows:

• We introduce the Safe Bayesian optimization with simulation-informed Gaussian process for the constraints (SI-SafeBO) algorithm, a variant of SafeOpt that enables its application with classification constraints. Although SafeOpt’s theoretical guarantees do not apply to classification constraints, we show that the proposed algorithm significantly reduces the number of constraint violations compared to directly learning classification constraints.

• We present an adaptation of SafeOpt that allows a GP classifier to directly learn true classification constraints, such as through Laplace approximation and expectation propagation. These methods are called safe Bayesian optimization with Laplace approximation (LA-SafeBO) and safe Bayesian optimization with expectation propagation (EP-SafeBO).

• We validate the algorithms’ performance using a computational model of a complex application, specifically the transport beamline of an ISOL system.

• We provide insights into the extent to which it is reasonable to use the available information from the constraints’ simulation model compared to directly learning the classification constraints, particularly considering the discrepancies between the simulation model and the true classification constraints.

The ISOL Technique

The production of radioisotopes and their delivery in the form of a radioactive ion beam (RIB) are carried out using the ISOL technique [45]. The process begins—at ISOL@MYRRHA—with the interaction between a proton beam, delivered by a linear accelerator (LINAC) with an energy of 100 MeV, and a target material. The target material is housed in the target and ion source assembly (TISA), typically heated up to 2300 °C and put to voltages up to 60 kV, as shown in Figure 2. The interaction between the proton beam and the target material leads to the formation of radioisotopes, illustrated as small dots. Additionally, unwanted isotopes (contaminants) are produced as by-products during irradiation. After evaporating from the target, these radioisotopes are ionized within the ion source (represented by dots with outer circles). The RIB is then formed by accelerating and extracting the ions due to the voltage difference between a grounded extraction electrode and the TISA (60 kV), represented by the shaded area in Figure 2. The RIB is transported from the extraction electrode (point A) to the RIB handling area. In this work, the analysis of the transport beamline focuses on the section up to the exit focal point of the dipole magnet (point J), as our primary interest lies in separating contaminants. Typically, experienced operators tune ISOL parameters for radioisotope extraction, but this process is time-consuming and may result in suboptimal performance. Consequently, optimization algorithms are crucial for accelerating the tuning process and achieving optimal performance.

At this point, it is pertinent to introduce the phase ellipse as the method for representing the RIB. The RIB comprises all the particles enclosed by the phase ellipse, defined in the six-dimensional phase space coordinates ($x,x^{\prime },y,y^{\prime },z,\Delta p/p$). Here, x and y represent the horizontal and vertical directions, whereas $x^{\prime }=\mathrm{d}x/\mathrm{d}z$ and $y^{\prime }=\mathrm{d}y/\mathrm{d}z$ denote the divergence or slope relative to the reference trajectory, z. $\Delta p/p$ indicates the relative momentum deviation from the ideal momentum p. The shape of the ellipse is characterized by the Twiss parameters ($\alpha ,\beta$, and $\gamma$) [46], with the emittance $ϵ$ defining the region occupied by particles within the beam in the phase space plane. It is important to note that the beam emittance exists in horizontal and vertical directions. A graphical representation of the phase ellipse is depicted in Figure 3, where $\mathrm{u}$ represents either the horizontal (x) or vertical (y) directions. Additionally, the beam envelope refers to the $2-\sigma$ boundary that encloses all particles within the beam, thereby defining the beam’s size as a function of the reference trajectory z.

2. Methods

2.1. Objective Function Definition

The optimization in this paper targets two goals: first, separating contaminants within the beam, and second, ensuring that the beam maintains opposite angles at the entry and exit focal points of the dipole magnet in the vertical direction. The latter goal helps keep the beam envelope consistent in size, preventing it from expanding and potentially reaching the physical limits of the magnet. The optimization problem is formulated to minimize the mathematical expression in Equation (1). The notation $(z=B\mid \mathbf{d})$ refers to the value of the Twiss parameters at the entrance focal point of the dipole magnet (point B from Figure 2), influenced by a set of quadrupole voltages $\mathbf{d}=\left\{d_1,d_2,...,d_8\right\}$. The reference Twiss parameters ${\beta }_x^{\mathrm{des}}(z=B)$, ${\alpha }_y^{\mathrm{des}}(z=B)$, and ${\beta }_y^{\mathrm{des}}(z=B)$ are values that must be reached to minimize the objective function; the first term from Equation (1) targets the separation of contaminants. The equation ${\sigma }_x=\sqrt{ϵ{\beta }_x}$ relates the horizontal beam size ${\sigma }_x$ and the Twiss parameter ${\beta }_x$. Therefore, the parameter ${\beta }_x(z=B\mid \mathbf{d})$ relates to the horizontal beam size at the dipole magnet’s entrance focal point. The choice of ${\beta }_x^{\mathrm{des}}(z=B)$ is based on a reference beam size that ensures contaminant separation up to a specific mass resolution or mass resolving power (MRP) [47]. Given the dipole magnet’s characteristics, if the beam attains the reference size, contaminants sufficiently deviate to be intercepted by a slit (point J in Figure 2). The remaining terms in Equation (1) ensure the angle condition is met at the dipole magnet’s entry and exit focal points in the vertical direction. Appendix A provides the derivation of the reference values ${\beta }_x^{\mathrm{des}}(z=B)$, ${\alpha }_y^{\mathrm{des}}(z=B)$, and ${\beta }_y^{\mathrm{des}}(z=B)$. The parameters $c_1,c_2,c_3$ represent weighting factors.

(1)$\begin{array}{ccc}\hfill f\left(\mathbf{d}\right)& =& c_1{[{\beta }_x(z=B\mid \mathbf{d})-{\beta }_x^{\mathrm{des}}(z=B)]}^2+c_2{[{\alpha }_y(z=B\mid \mathbf{d})-{\alpha }_y^{\mathrm{des}}(z=B)]}^2\hfill \\ & +& c_3{[{\beta }_y(z=B\mid \mathbf{d})-{\beta }_y^{\mathrm{des}}(z=B)]}^2.\hfill \end{array}$

Ensuring that the beam envelope stays within the physical limits of the beamline during transport is crucial for its successful delivery to the intended destination. Therefore, constraints are formulated to meet this requirement. The transport beamline constraints from location A to G are defined as $g^i\left(\mathbf{d}\right)=\max \left[{\sigma }_{x,y}(z\mid \mathbf{d})\right]-R_{A{G}_{x,y}}\left(z\right)$ for all $i=1,...,n\mathrm{and}A\le z\le G$. There are elements with different aperture size along the transport beamline, and the constraints are governed by the element with the smallest aperture ($R_{A{G}_{x,y}}\left(z\right)$). The parameter tuning process of the ISOL system is therefore framed as a multidimensional constrained optimization problem, formulated as

(2)$\underset{\mathbf{d}}{\min }f\left(\mathbf{d}\right)\mathrm{subject}\mathrm{to}{g}^i\left(\mathbf{d}\right)<0\mathrm{for}i=1,...,4.$

The ISOL system’s transport beamline instrumentation allows for measuring the beam size at specific locations ($z=A,B,J$). Consequently, the term $\max \left[{\sigma }_{x,y}(z\mid \mathbf{d})\right]$ in the expression $g^i\mathrm{for}i=1,...,4$ cannot be directly computed. Instead, constraint analysis is simplified to a classification task, determining whether the beam successfully reaches its destination. If the beam reaches its destination, it is labeled safe and represented numerically by ${\mathrm{y}}^i=-1$. Conversely, if the beam collides with a beamline element, the constraint is labeled unsafe and represented by ${\mathrm{y}}^i=+1$.

2.2. Gaussian Process

The Gaussian process (GP) is a widely used approach for defining models of nonlinear functions [48]. A significant advantage of GP is its ability to provide a probabilistic characterization of function predictions at each input location, offering both the predicted function value and the associated uncertainty. This section describes the general procedure for constructing a GP model to predict an unknown mathematical function $\mathrm{y}=f\left(\mathbf{d}\right)+\nu$, where $\nu$ represents zero-mean Gaussian-distributed noise with variance ${\sigma }_{\mathbf{y}}^2$. In this context, $\mathbf{d}=(d_1,...,d_r)$ denotes an input location. $\mathcal{F}=\left\{\mathbf{D},\mathbf{y}\right\}$ represents a set of N observations from the function, where $\mathbf{y}=\left\{{\mathrm{y}}_1;...;{\mathrm{y}}_N\right\}\in {\mathbb{R}}^N,\mathrm{and}\mathbf{D}=\left\{{\mathbf{d}}_1;...;{\mathbf{d}}_N\right\}\in {\mathbb{R}}^{N\times r}$. A GP can be understood as a generator of numbers using a multivariate normal distribution $\mathbf{y}\sim \mathcal{N}(\mathit{\mu },\mathsf{\Sigma })$ with mean $\mathit{\mu }=\left({\mu }_1,...,{\mu }_N\right)$ and covariance matrix $\mathsf{\Sigma }=\mathbf{K}(\mathbf{D},\mathbf{D})+{\sigma }_{\mathbf{y}}^2\mathbf{I}$ to generate a vector of N numbers, $\mathbf{y}=\left({\mathrm{y}}_1;...;{\mathrm{y}}_N\right)$. A GP constructs the matrix $\mathbf{K}(\mathbf{D},\mathbf{D})$ by introducing covariance between successive input locations $\left(\mathbf{d},{\mathbf{d}}^{\prime }\right)$ via a kernel function $k(\mathbf{d},{\mathbf{d}}^{\prime })$, such as the squared exponential kernel [49]:

(3)$k(\mathbf{d},{\mathbf{d}}^{\prime })={\sigma }_f^2\exp \left[-0.5\left(\mathbf{d}-{\mathbf{d}}^{\prime }\right)\mathsf{\Lambda }{\left(\mathbf{d}-{\mathbf{d}}^{\prime }\right)}^T\right].$

The kernel function parameters, along with ${\sigma }_{\mathbf{y}}^2$, are referred to as the hyperparameters, denoted as $\mathit{\lambda }=\left\{\mathbf{l},{\sigma }_f,{\sigma }_{\mathbf{y}}^2\right\}$ when using Equation (3), where $\mathsf{\Lambda }=\mathrm{diag}{\left(\mathbf{l}\right)}^{-2}\in {\mathbb{R}}^{r\times r}$ and $\mathbf{l}={\mathrm{l}}_1,...,{\mathrm{l}}_r$ allowing each variable to have an associated length scale $\mathrm{l}$. The numbers generated from each sample of the multivariate normal distribution $\mathbf{y}\sim \mathcal{N}(\mathit{\mu },\mathsf{\Sigma })$ represent a function realization. Using the observed input locations $\mathbf{D}$ and a set of hyperparameters $\mathit{\lambda }$, a sample from the multivariate distribution $\mathcal{N}(\mathit{\mu }\left(\mathbf{D}\right)=\mathbf{0},\mathsf{\Sigma }=\mathbf{K}(\mathbf{D},\mathbf{D}|\mathit{\lambda })+{\sigma }_{\mathbf{y}}^2\mathbf{I})$ generates a function realization that estimates the unknown mathematical function $f\left(\mathbf{D}\right)$. The second step in the GP modeling procedure, known as hyperparameter selection, involves selecting hyperparameters that produce a realization closely matching the function $f\left(\mathbf{d}\right)$. Optimal hyperparameters for GP are obtained by ${\mathit{\lambda }}^{\mathrm{opt}}={\mathrm{argmax}}_{\mathit{\lambda }}logp[\mathbf{y}\mid \mathbf{D},\mathit{\lambda }]$, where the expression $p[\mathbf{y}\mid \mathbf{D},\mathit{\lambda }]$ represents the probability of obtaining the observed function values $\mathbf{y}$, given the input locations $\mathbf{D}$ and the hyperparameters $\mathit{\lambda }$. Therefore, sampling the multivariate normal distribution $\mathcal{N}(\mathit{\mu }\left(\mathbf{D}\right)=\mathbf{0},\mathbf{K}(\mathbf{D},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}})+{\sigma }_{\mathbf{y}}^2\mathbf{I})$ with the optimal hyperparameters ${\mathit{\lambda }}^{\mathrm{opt}}$ generates function realizations that best capture the observed function values $\mathbf{y}$.

The final step in GP modeling involves conditioning the multivariate normal distribution on the observations $\mathbf{y}$ to produce a realization that predicts the function ${\mathbf{y}}_{*}\sim f\left({\mathbf{D}}_{*}\right)$ at query points ${\mathbf{D}}_{*}=\left\{{\mathbf{d}}_{*1};...;{\mathbf{d}}_{*M}\right\}\in {\mathbb{R}}^{M\times r}$, which may differ from the observed input locations $\mathbf{D}$. The latter is possible because the GP formulation assumes that the joint distribution of ${\left[\begin{array}{cc}\mathbf{y}& {\mathbf{y}}_{*}\end{array}\right]}^T$ can be approximated by the multivariate normal distribution:

(4)$\left[\begin{array}{c}\mathbf{y}\\ {\mathbf{y}}_{*}\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\end{array}\right]\left[\begin{array}{cc}\mathbf{K}(\mathbf{D},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}})+{\sigma }_{\mathbf{y}}^2\mathbf{I}& \mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}|{\mathit{\lambda }}^{\mathrm{opt}})\\ \mathbf{K}({\mathbf{D}}_{*},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}})& \mathbf{K}({\mathbf{D}}_{*},{\mathbf{D}}_{*}|{\mathit{\lambda }}^{\mathrm{opt}})\end{array}\right]\right).$

The posterior distribution of ${\mathbf{y}}_{*}$ is then obtained, as expressed in Equation (5), characterized by its mean $\mathit{\mu }\left({\mathbf{D}}_{*}\right)$ and variance $\mathsf{\Sigma }\left({\mathbf{D}}_{*}\right)$. Thus, after selecting the hyperparameters as described earlier, the conditioning step determines the prediction’s mean. Moreover, the posterior distribution at each query point ${\mathbf{d}}_{*j}$ is approximated as a univariate normal distribution. The mean of this distribution equals the jth element of the vector $\mathit{\mu }\left({\mathbf{D}}_{*}\right)$, and its variance corresponds to the jth element on the diagonal of the matrix ${\sigma }^2\left({\mathbf{D}}_{*}\right)=\mathsf{\Sigma }{\left({\mathbf{D}}_{*}\right)}_{i,j},\forall i=j\in {\mathbb{N}}_{[1,M]}$. For clarity, we denote the posterior distribution as ${\mathcal{GP}}^{\mathit{\mu }}=\mathit{\mu }\left({\mathbf{D}}_{*}\right)$ for the mean and ${\mathcal{GP}}^{\sigma }=\sqrt{\mathsf{\Sigma }\left({\mathbf{D}}_{*}\right)}$ for the standard deviation.

(5)$\begin{array}{ccc}\hfill f\left({\mathbf{D}}_{*}\right)& |& \mathbf{D},\mathbf{y},{\mathbf{D}}_{*}\sim \mathcal{N}\left(\mathit{\mu }\left({\mathbf{D}}_{*}\right),\mathsf{\Sigma }\left({\mathbf{D}}_{*}\right)\right),\hfill \\ \hfill \mathit{\mu }\left({\mathbf{D}}_{*}\right)& =& \mathbf{K}({\mathbf{D}}_{*},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}}){\left[\mathbf{K}(\mathbf{D},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}})+{\sigma }_{\mathbf{y}}^2\mathbf{I}\right]}^{-1}\mathbf{y},\hfill \\ \hfill \mathsf{\Sigma }\left({\mathbf{D}}_{*}\right)& =& \mathbf{K}({\mathbf{D}}_{*},{\mathbf{D}}_{*}|{\mathit{\lambda }}^{\mathrm{opt}})-\mathbf{K}({\mathbf{D}}_{*},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}}){\left[\mathbf{K}(\mathbf{D},\mathbf{D}|{\mathit{\lambda }}^{\mathrm{opt}})+{\sigma }_{\mathbf{y}}^2\mathbf{I}\right]}^{-1}\hfill \\ & & \times \mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}|{\mathit{\lambda }}^{\mathrm{opt}}).\hfill \end{array}$

2.3. Safe Bayesian Optimization

Sui et al. [27] introduced safe Bayesian optimization (SafeOpt), a formulation that employs Bayesian optimization (BO) to find an optimal solution while ensuring constraint satisfaction throughout the optimization process. This paper extends SafeOpt to handle classification constraints, provided a simulation model of those constraints is available. The proposed formulation addresses the absence of the Lipschitz constant L in the construction of the safe set, as proposed by Berkenkamp et al. [30], and incorporates the line search approach from line Bayesian optimization [14]. The resulting algorithm is referred to as SafeBO. In this context, the superscript zero denotes the objective function, while superscripts from one to n refer to the constraints. The lower and upper confidence bounds are denoted as $l_t^i={\mathcal{GP}}_t^{{\mu }_i}-\kappa {\mathcal{GP}}_t^{{\sigma }_i}$ and $u_t^i={\mathcal{GP}}_t^{{\mu }_i}+\kappa {\mathcal{GP}}_t^{{\sigma }_i}\mathrm{for}\mathrm{all}i=0,...,n$, respectively. In this work, we adopt the notation $\kappa$ for the scaling factor [50] instead of the standard notation $\beta$ used in the Gaussian process literature to avoid confusion with the Twiss parameter from Section 1. The width of the bounds is defined as $w_t^i=u_t^i-l_t^i$. The first step involves constructing the safe set $\mathcal{S}$, as illustrated in the one-dimensional example with $r=1$ and one constraint $i=0,1$ shown in Figure 4a. The mathematical expressions in this section are given for the general case $r\ge 1$. The algorithm begins with an initial safe set ${\mathcal{S}}_0$, represented by the dot. SafeBO then proceeds to build the safe set using input locations $\mathbf{d}\in {\mathcal{A}}_t$ where the upper confidence bounds of the GP prediction for the constraints fall below a negative safety threshold $\tau$,

(6)${\mathcal{S}}_t\leftarrow \bigcap _{i\in 1,...,n}\bigcup _{\mathbf{d}\in {\mathcal{S}}_{t-1}}\left\{\mathbf{d}\in {\mathcal{A}}_t\mid u_{t-1}^i\left(\mathbf{d}\right)\le \tau \right\}.$

The remaining locations ${\mathbf{d}}^{\prime }\in {\mathcal{A}}_t\setminus {\mathcal{S}}_t$ are considered unsafe. The safe set search is restricted to a one-dimensional subspace ${\mathcal{A}}_t$ to avoid intractability in high-dimensional spaces, as suggested by Kirschner et al. [28]. This subspace contains the current best solution (${\mathbf{d}}_{t-1}^{+}$) and defines a line search in a specific direction b. SafeBO randomly selects locations within a ball $\mathcal{B}({\mathbf{d}}_{t-1}^{+},\eta )$ with radius $\eta$, choosing the direction based on where the objective function indicates the most significant potential decrease. The second step involves constructing the set of possible minimizers $\mathcal{M}$, as shown in Figure 4b. This set includes only the locations from the safe set $\mathcal{S}$ where the lower bound of the objective function $l_{t-1}^0$ is less than the location that attains the lowest upper bound

(7)${\mathcal{M}}_t\leftarrow \left\{\mathbf{d}\in {\mathcal{S}}_t\mid l_{t-1}^0\left(\mathbf{d}\right)\le \underset{{\mathbf{d}}^{\prime }\in {\mathcal{S}}_t}{\min }{u}_{t-1}^0\left({\mathbf{d}}^{\prime }\right)\right\}.$

SafeBO is rooted in Bayesian optimization and retains the core principle of balancing exploitation and exploration. Exploitation, leveraging the model’s current knowledge, is focused on the set of minimizers $\mathcal{M}$. At the same time, exploration is achieved by constructing a set of potential expanders $\mathcal{G}$, which further expand the safe set $\mathcal{S}$, allowing for cautious exploration of the design space. Figure 4c illustrates the method for constructing the set of expanders, which includes locations from $\mathcal{S}$ that satisfy the optimistic indicator function [30].

(8)${\mathcal{G}}_t\leftarrow \left\{\mathbf{d}\in {\mathcal{S}}_t,\mid e_t^i\left(\mathbf{d}\right)\le 0\right\}.$

To understand the expression in the optimistic indicator function, we refer to the star in the lower part of Figure 4c. This star represents an artificial observation consisting of a safe location $\mathbf{d}$ and its lower confidence bound $l_{t-1}^i$. The idea is to update the GP models of the constraints by incorporating this artificial location along with past observations, forming ${\mathcal{GP}}_t^{a,i}\mid {\mathcal{C}}_{t-1}\cup \left(\mathbf{d},l_{t-1}^i\left(\mathbf{d}\right)\right)$. Then, the upper confidence bounds are computed at the unsafe locations using the updated GP predictions for the constraints, defined as $u_{t,\left(\mathbf{d},l_{t-1}^i\left(\mathbf{d}\right)\right)}^i\left({\mathbf{d}}^{\prime }\right)={\mathcal{GP}}_t^{a,{\mu }_i}\left({\mathbf{d}}^{\prime }\right)+\kappa {\mathcal{GP}}_t^{a,{\sigma }_i}\left({\mathbf{d}}^{\prime }\right)$, so that

(9)$e_t^i\left(\mathbf{d}\right)=\left\{{\mathbf{d}}^{\prime }\in {\mathcal{A}}_t\setminus {\mathcal{S}}_t\mid u_{t,\left(\mathbf{d},l_{t-1}^i\left(\mathbf{d}\right)\right)}^i\left({\mathbf{d}}^{\prime }\right)\le \tau \right\}.$

The safe location $\mathbf{d}$ is considered a potential expander if the predicted upper bound at unsafe locations falls below the threshold $\tau$. Including this safe location in the next iteration reduces uncertainty near the unsafe locations, helping to identify the feasibility limits of the design space. SafeBO’s final decision for selecting the best location ${\mathbf{d}}_t^{*}$ is made by balancing exploration and exploitation, choosing the location from $\mathcal{G}\cup \mathcal{M}$ that exhibits the widest uncertainty bounds in the GP predictions for the objective function and constraints.

(10)${\mathbf{d}}_t^{*}\leftarrow \underset{\mathbf{d}\in {\mathcal{G}}_t\cup {\mathcal{M}}_t}{\mathrm{argmax}}\underset{i\in 1,...,n}{\max }{w}_{t-1}^i\left(\mathbf{d}\right).$

In the illustrative example from Figure 4d, the proposed best location achieves the widest uncertainty bound for the objective function prediction, as no location shows wider uncertainty bounds for the constraint prediction. In this paper, we use normalized values for the objective function f and the input locations $\mathbf{d}$ within the range [0–1], while the constraints $g^i$ are normalized within the range [−1–0]. This normalization ensures a fair comparison of the GP model’s confidence bounds for the objective function and constraints. Kirschner et al. [14] emphasize the importance of normalizing constraints to negative values, as this allows for identifying unsafe input parameters in the absence of data. Moreover, normalizing the constraint results in a negative threshold $\tau$ that defines the safety condition. However, adjusting this threshold can introduce a safety margin, which may slow convergence [14].

2.3.1. Safe Bayesian Optimization for Classification Constraints

Classification constraints can be modeled using techniques such as Laplace approximation and expectation propagation [41]. Both approaches use a latent function $h\left(\mathbf{d}\right)$, an intermediate valuable function for model formulation. The key idea behind Gaussian process approximation is to fit a GP prior over the latent function and then map it to the [0–1] domain using the logistic or probit function, as illustrated in Figure 5. This mapping results in the prior distribution on $p(\mathrm{y}=+1|\mathbf{d})=\zeta (h\left(\mathbf{d}\right))$, where $\zeta (.)$ denotes the logistic or probit function. Predictions require the posterior distribution, which involves first computing the latent variable distribution for the query locations ${\mathbf{D}}_{*}$ and subsequently using this distribution to compute the probabilistic prediction:

(11)$\begin{array}{ccc}\hfill p\left(h\left({\mathbf{D}}_{*}\right)\right|\mathbf{D},\mathbf{y},{\mathbf{D}}_{*})& =& \int p\left(h\left({\mathbf{D}}_{*}\right)\right|\mathbf{D},{\mathbf{D}}_{*},h\left(\mathbf{D}\right))p\left(h\left(\mathbf{D}\right)\right|\mathbf{D},\mathbf{y})dh\left(\mathbf{D}\right),\hfill \\ \hfill p({\mathbf{y}}_{*}=+\mathbf{1}|\mathbf{D},\mathbf{y},{\mathbf{D}}_{*})& =& \int \zeta \left(h\left({\mathbf{D}}_{*}\right)\right)p\left(h\left({\mathbf{D}}_{*}\right)\right|\mathbf{D},\mathbf{y},{\mathbf{D}}_{*}).\hfill \end{array}$

The inconvenience with the classification case is that the integral in Equation (11) becomes analytically intractable because the nonlinear mapping introduced by $\zeta (.)$ makes the joint distribution $p\left(\mathbf{y}\right|h\left(\mathbf{D}\right))$ non-Gaussian. Therefore, analytical approximations of these integrals are obtained through the Laplace approximation or expectation propagation. The Laplace approximation uses Newton’s method to generate a second-order Taylor expansion around the mode $\mathbf{m}$. The posterior distribution over the latent function and the probabilistic prediction are expressed as

(12)$\begin{array}{ccc}\hfill \mathit{\mu }\left({\mathbf{D}}_{*}\right)& =& \mathbf{K}(\mathbf{D},{\mathbf{D}}_{*})\nabla \log p\left(\mathbf{y}\right|h\left(\mathbf{D}\right)),\hfill \\ \hfill \mathsf{\Sigma }\left({\mathbf{D}}_{*}\right)& =& \mathbf{K}({\mathbf{D}}_{*},{\mathbf{D}}_{*})-\mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}){\left[\mathbf{K}(\mathbf{D},\mathbf{D})+{[-\nabla \nabla \log p\left(\mathbf{y}\right|h\left(\mathbf{D}\right))]}^{-1}\right]}^{-1}\hfill \\ & \times & \mathbf{K}({\mathbf{D}}_{*},\mathbf{D}),\hfill \\ \hfill p({\mathbf{y}}_{*}=+\mathbf{1}|\mathbf{D},\mathbf{y},{\mathbf{D}}_{*})& =& \int \zeta \left(h\left({\mathbf{D}}_{*}\right)\right)\mathcal{N}(\mathit{\mu }\left({\mathbf{D}}_{*}\right),\mathsf{\Sigma }\left({\mathbf{D}}_{*}\right))dh\left({\mathbf{D}}_{*}\right).\hfill \end{array}$

For expectation propagation, the likelihood is approximated by a local approximation that uses an un-normalized Gaussian on the latent function with site parameters $\tilde{\mathit{\mu }},\tilde{\mathsf{\Sigma }}$. In this case, the posterior distribution and the probabilistic prediction are

(13)$\begin{array}{ccc}\hfill \mathit{\mu }\left({\mathbf{D}}_{*}\right)& =& \mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}){[\mathbf{K}(\mathbf{D},\mathbf{D})+\tilde{\mathsf{\Sigma }}]}^{-1}\tilde{\mathit{\mu }},\hfill \\ \hfill \mathsf{\Sigma }\left({\mathbf{D}}_{*}\right)& =& \mathbf{K}({\mathbf{D}}_{*},{\mathbf{D}}_{*})-\mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}){\left[\mathbf{K}(\mathbf{D},\mathbf{D})+\tilde{\mathsf{\Sigma }}\right]}^{-1}\mathbf{K}({\mathbf{D}}_{*},\mathbf{D}),\hfill \\ \hfill p({\mathbf{y}}_{*}=+\mathbf{1}|\mathbf{D},\mathbf{y},{\mathbf{D}}_{*})& =& \mathsf{\Phi }\left({\displaystyle \frac{\mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}){[\mathbf{K}(\mathbf{D},\mathbf{D})+\tilde{\mathsf{\Sigma }}]}^{-1}\tilde{\mathit{\mu }}}{\sqrt{1+\mathbf{K}({\mathbf{D}}_{*},{\mathbf{D}}_{*})-\mathbf{K}(\mathbf{D},{\mathbf{D}}_{*}){\left[\mathbf{K}(\mathbf{D},\mathbf{D})+\tilde{\mathsf{\Sigma }}\right]}^{-1}\mathbf{K}({\mathbf{D}}_{*},\mathbf{D})}}}\right).\hfill \end{array}$

Further details about the mathematical derivation and implementation of the Laplace approximation and expectation propagation can be found in Nickisch and Rasmussen [24], Rasmussen and Williams [41]. The application of SafeBO with classification constraints necessitates a modification in the construction of the set of expanders $\mathcal{G}$. Instead of using the lower confidence bound $(\mathbf{d},l_{t-1}^i\left(\mathbf{d}\right))$, the artificial location incorporates the probabilistic classification for the constraints $(\mathbf{d},{\mathrm{y}}_t^{a,i})$. This adjustment is necessary because using the lower confidence bound at the artificial location might include many infeasible input locations in the safe set, given that the GP approximation for classification typically results in high uncertainty bounds. The probabilistic classification is based on the probability of $\mathbf{d}$ belonging to the safe or unsafe class, as derived in Equation (12) for the Laplace approximation and Equation (13) for the expectation propagation. The artificial input location is labeled safe if $p({\mathrm{y}}_t^{a,i}=+1\mid \mathbf{d})$ remains below $50\%$; otherwise, it is labeled unsafe. In the Laplace approximation and expectation propagation, the value −1 denotes a safe location, while 1 indicates an unsafe location, consistent with the GPML MATLAB toolbox requirements [41].

Figure 6 presents an illustrative example highlighting the application of SafeBO with Laplace approximation for constraints (LA-SafeBO) to minimize an objective function, where the objective function is omitted for clarity. Furthermore, the predicted values for the constraints are normalized to [−1–0] after classification to enable a fair comparison among the confidence bounds, with −1 indicating safe and 0 indicating unsafe. In Figure 6a, the dots represent the constraint observations selected by the algorithms after 20 iterations. The solid line depicts the constraint function $g^1\left(d\right)$ (ground truth). The dashed line and shaded area represent the GP approximation for the constraint with its mean ${\mathcal{GP}}^{{\mu }_1}$ and standard deviation ${\mathcal{GP}}^{{\sigma }_1}$, following the notation from the previous section, where superscript 1 refers to the constraint. The dash-dotted line shows the probabilistic classification of the constraint ${\widehat{g}}^1\left(d\right)$.

Notably, the safe set $\mathcal{S}$ includes an infeasible region of the design space. This example illustrates the challenge of learning discontinuous functions as the GP approximation struggles to predict an unexpected change in the function. The probabilistic classification shows that the entire design space is classified as safe before a constraint violation occurs. This is logical since only safe locations are included in the training data. Thus, as safeBO uses the upper bound of the GP approximation to build the safe set $\mathcal{S}$, it may include infeasible input locations, resulting in a location $d^{*}$ that violates the constraint. Additionally, Figure 6b shows the final model of the constraint after 30 iterations. The Laplace approximation learns the feasibility limits of the constraint with relatively good accuracy after violating it on four occasions during the learning process. However, despite the multiple observations, the GP approximation still exhibits wide uncertainty bounds. Finally, this highlights that learning constraints may require observations in unsafe regions, which is suboptimal for safety applications.

2.3.2. Safe Bayesian Optimization with Simulation-Informed Gaussian Process for the Constraints

This section introduces safe Bayesian optimization with simulation-informed GP for the constraints (SI-SafeBO), a variation of SafeBO designed to reduce observations in unsafe regions when dealing with classification constraints. Since SI-SafeBO builds on SafeBO, it begins by selecting the input location most likely to minimize the objective function while satisfying the constraints (${\mathbf{d}}_t^{*}$), as described in Section 2.3. The objective function and constraints are then evaluated at this location. However, SI-SafeBO proposes a constraint function ($g^{p,i}$) that combines observations from both the true classification constraint and the simulation constraint. Rather than assigning a numerical value of −1 for safe locations, SI-SafeBO utilizes continuous observations from the simulation constraint, such as $d=[0.06–0.43]\cup [0.57–0.86]$ in the one-dimensional example from Figure 7a. Conversely, when a constraint violation occurs (${\mathrm{y}}_t^i$ takes the value 0), SI-SafeBO uses the ground truth constraint value ($g^i$), represented by regions with a zero value in Figure 7a. As a result, the constraint for SI-SafeBO ($g^{p,i}$) becomes a piece-wise function, depicted by the solid line in Figure 7a.

In scenarios where no constraint violations occur, the GP of the constraint models the continuous segment of $g^{p,i}$ accurately. In Figure 7a, the dashed line and shadow represent the GP’s prediction for the constraint after two iterations. The triangle represents the SafeBO suggested location $d^{*}$ for the next iteration, and the dotted line marks the safety threshold. It can be seen that the GP is learning the continuous segment of the piece-wise function. Figure 7b,c depict the GP of the constraint after 13 and 30 iterations, respectively. These figures show that the uncertainty bounds for the simulation-informed GP (SI-GP) are narrower than those for the Laplace approximation. This effect is attributed to the continuous segment of the piece-wise constraint because the observations follow a continuous function. Furthermore, if a constraint violation occurs, the simulation-informed GP can still approximate the piece-wise function, as depicted in Figure 7d. Since the GP uses continuous observations from the simulation constraint, it attempts to fit a continuous function to the piece-wise function, allowing it to accurately learn the limits of the feasible design space. In the example shown in Figure 7d, it is clear that the updated upper bound of the SI-GP is above the safety threshold at locations where the infeasible design space begins. In summary, the SI-GP approach replaces the true classification constraint with a piece-wise function, where the continuous part comes from a simulation constraint. This approach enables the SI-GP to approximate the piece-wise function and apply SafeBO with classification constraints.

However, discrepancies between the simulation and true classification constraints may arise, potentially leading SI-SafeBO to suggest a location that violates the simulation constraint but remains safe in the true constraint. This situation is illustrated in Figure 8a. The suggested location $d^{*}$, represented by the triangle, is in an infeasible region of the piece-wise function. Nevertheless, it is evident that at location $d^{*}$, the true constraint function $g^1$ is feasible. To address this issue, SI-GP relies on the GP’s predicted mean value $\left({\mathcal{GP}}^{{\mu }_i}\right)$ of the constraint at that location, represented by the star. This approach helps circumvent the problem of learning incorrect feasible region limits. In Figure 8b, the SI-SafeBO algorithm learned the true feasible limits of the constraint despite the discrepancies between the simulation and true constraints. In this example, only one constraint violation was required to learn the true feasibility limits, compared to the four violations needed with the Laplace approximation.

It is worth noting that constraint violations may still occur due to discrepancies between the simulation model and the ground truth system. However, the infeasible region is likely to be avoided in future observations as the constraint’s GP model updates and learns the boundaries of the infeasible region. The advantages of the proposed approach in reducing the number of constraint violations compared to directly learning the true classification constraints using the Laplace approximation or expectation propagation are emphasized in the results section. This method applies only to processes where a simulation constraint exists. Lastly, Algorithm 1 provides the pseudo-code of the SI-SafeBO algorithm.

Algorithm 1 Safe Bayesian optimization with simulation-informed GP for the constraints

Require:  Initial safe set ${\mathcal{S}}_0,{\mathcal{F}}_0,{\mathcal{C}}_0^i$, best initial location ${\mathbf{d}}_0^{+}$,   . objective function and constraint GP priors $\mathcal{N}({\mathcal{GP}}_0^{{\mu }_i},{\mathcal{GP}}_0^{{\mathit{\sigma }}_i})$   . $\forall i\in 0,\dots ,n$ 1:  . for $t=1,2,3,\dots$do 2:.    ${\mathbf{d}}_t^{*}\leftarrow \mathrm{SafeBO}({\mathbf{d}}_{t-1}^{+},{\mathcal{S}}_{t-1},\mathcal{N}({\mathcal{GP}}_{t-1}^{{\mu }_i},{\mathcal{GP}}_{t-1}^{{\mathit{\sigma }}_i}))$               ▹ from Section 2.3 3:.      ${\mathrm{y}}_t^0\leftarrow f_t\left({\mathbf{d}}_t^{*}\right)+{\nu }^0$ 4:.      ${\mathrm{y}}_t^i\leftarrow g_t^i\left({\mathbf{d}}_t^{*}\right)+{\nu }^i\forall i\in 1,\dots ,n$ 5:.      if ${\mathrm{y}}_t^i<0$ then 6:.        ${\mathrm{y}}_t^i\leftarrow g_t^{p,i}\left({\mathbf{d}}_t^{*}\right)$ 7:.        if ${\mathrm{y}}_t^i\ge 0$ then 8:.         ${\mathrm{y}}_t^i\leftarrow {\mathcal{GP}}_{t-1}^{{\mu }_i}\left({\mathbf{d}}_t^{*}\right)$ 9:.        end if 10:.      end if 11:.      ${\mathcal{F}}_t\leftarrow {\mathcal{F}}_{t-1}\cup ({\mathbf{d}}_t^{*},{\mathrm{y}}_t^0)$ 12:.      ${\mathcal{C}}_t^i\leftarrow {\mathcal{C}}_{t-1}^i\cup ({\mathbf{d}}_t^{*},{\mathrm{y}}_t^i)$ 13:.      $\mathcal{N}({\mathcal{GP}}_t^{{\mu }_0},{\mathcal{GP}}_t^{{\mathit{\sigma }}_0})\leftarrow \mathcal{N}({\mathcal{GP}}_{t-1}^{{\mu }_0},{\mathcal{GP}}_{t-1}^{{\mathit{\sigma }}_0})\mid {\mathcal{F}}_t$ 14:.      $\mathcal{N}({\mathcal{GP}}_t^{{\mu }_i},{\mathcal{GP}}_t^{{\mathit{\sigma }}_i})\leftarrow \mathcal{N}({\mathcal{GP}}_{t-1}^{{\mu }_i},{\mathcal{GP}}_{t-1}^{{\mathit{\sigma }}_i})\mid {\mathcal{C}}_t^i$ 15:.      if ${\mathrm{y}}_t^0<f\left({\mathbf{d}}_{t-1}^{+}\right)$ then 16:.      ${\mathbf{d}}_t^{+}\leftarrow {\mathbf{d}}_t^{*}$ 17:.      end if 18:  . end for 19:  . ${\mathbf{d}}^{\mathrm{opt}}\leftarrow {\mathbf{d}}^{+}$ 20:  . return ${\mathbf{d}}^{\mathrm{opt}}$

3. Results and Discussion

This section presents the results of solving the optimization problem from Equation (2) for a computational model of the ISOL@MYRRHA facility’s transport beamline (ground truth). This model classifies constraints as safe (${\mathrm{y}}^i=-1$) or unsafe (${\mathrm{y}}^i=+1$). An additional simulation model provides continuous constraint values used to build the constraints $g^{p,i}$ employed by the simulation-informed GP. The computational model and optimization algorithms were implemented in MATLAB R2021b.

We compare SafeBO with simulation-informed GP for the constraints (SI-SafeBO), SafeBO with Laplace approximation (LA-SafeBO), and SafeBO with expectation propagation (EP-SafeBO). A GP models the real-valued objective function for each method, as depicted in Figure 1. Consequently, the method’s name reflects the approach to approximate the classification constraints. The algorithms are evaluated based on two performance criteria: convergence speed and the number of constraint violations. The optimization starts from 10 different voltage combinations that satisfy all constraints, simulating an operator’s tuning to a suboptimal solution. The optimization is run five times for each starting voltage combination, resulting in 50 runs. Each run is limited to 80 function evaluations, sufficient to reach the objective function threshold ($f^{\mathrm{trh}}$), which determines convergence. Notice that the computational effort required to compute the objective function is independent of the optimization algorithms. Therefore, the convergence speed is determined by the number of objective function evaluations the algorithm needs to reach the objective threshold. Measurement device noise is considered by adding Gaussian-distributed noise with ${\sigma }_{\mathbf{y}}^2=0.1$ to the objective function. Additionally, critical parameters for SafeBO must be determined at the outset of the optimization.

These parameters include $\kappa$, $\eta$, and $\tau$ from Section 2.3. In SafeBO, the construction of the safe set $\mathcal{S}$ depends on $\kappa$ so that a larger $\kappa$ leads to more conservative upper and lower bounds. In contrast, a smaller value can result in constraint violations, especially early on when uncertainty is high. The parameter $\eta$ represents the extent of the line search and is related to constraint violations because a larger value allows for selecting locations with higher uncertainty. $\tau$ acts as a safety factor by lowering the safety threshold by subtracting a percentage of $\tau$, which may reduce constraint violations but lead to a more conservative algorithm and slower convergence. In this paper, we studied the effect of these parameters on the algorithm’s convergence and the number of constraint violations. Based on the results in Appendix B, the selected values for the parameters are $\kappa =5.5$, $\eta =0.1$, and $\tau =0\%$. A critical objective function value defines the objective threshold ($f^{\mathrm{trh}}$), where the probability of achieving a mass resolving power (MRP) of 300 is approximately $85\%$; this specific value is justified in Appendix A. It is worth noting that the current of the dipole magnet can be adjusted, and depending on the set current, a specific mass passes through the reference trajectory. In the experiments of this section, we assume that the current of the dipole magnet is fixed. Constraint violations occur when the beam envelopes reach the physical limits of the transport beamline elements. Figure 9a–c presents the convergence results. The horizontal axis represents the number of objective function evaluations, whereas the vertical axis represents the normalized objective function f from Equation (1). The thick lines indicate the median of f across 50 independent optimizations, and the shaded areas represent the 10th and 90th percentiles. The dotted line indicates the objective threshold.

All three algorithms require approximately 25 objective function evaluations to reach the objective threshold, demonstrating similar convergence speeds. However, the median of the objective function value for LA-SafeBO and EP-SafeBO exhibits steeper slopes towards the objective threshold than SI-SafeBO because the GP classifier predicts that the design space is feasible until it learns the infeasible regions, allowing more exploration. In contrast, the simulation model used in SI-SafeBO allows for a cautious exploration of the design space. Considering that all algorithms use SafeBO to select the proposed best location ${\mathbf{d}}_t^{*}$ at each iteration, they were expected to achieve similar convergence. Figure 9d presents a box plot showing the variation in constraint violations for each algorithm across all 50 optimizations. SI-SafeBO achieved a median of two violations, whereas LA-SafeBO recorded a median of 10 violations, and EP-SafeBO had a median of $11.5$. The higher number of constraint violations observed with LA-SafeBO and EP-SafeBO result from directly learning the classification constraints. At the beginning of the optimization, the probabilistic prediction tends to classify much of the design space as feasible, allowing for exploration until the algorithm learns the feasibility limits. Additionally, the GP approximation of the constraints ($({\mathcal{GP}}^{{\mathit{\mu }}_i},{\mathcal{GP}}^{{\sigma }_i})\forall i=1,...,n$) for LA-SafeBO and EP-SafeBO generates high uncertainty in the predictions, necessitating more information from the infeasible design space to accurately learn the boundaries of the feasible design space. Furthermore, in $16\%$ of the runs, SI-SafeBO completes the optimization routine without violating any constraints, compared to $2\%$ for LA-SafeBO and $4\%$ EP-SafeBO, showing another advantage in using the proposed method.

Considering that the three algorithms achieve relatively similar convergence speeds, they require the same effort in terms of objective function evaluations. However, another relevant aspect is the time SI-SafeBO, LA-SafeBO, and EP-SafeBO take to suggest the best location ${\mathbf{d}}_t^{*}$ at each iteration. These times were measured for one run out of the 50 runs conducted in the experiment of this section, where the mean and maximum values over the 80 iterations were computed. Table 1 summarizes the results. In this regard, SI-SafeBO shows a slight advantage over LA-SafeBO, while EP-SafeBO requires the most time to suggest a solution. Based on the results in this section, SI-SafeBO is the recommended algorithm for optimizing applications with classification constraints when a simulation model of the constraints is available. If a simulation model is unavailable, LA-SafeBO is suggested, as it achieves similar convergence and constraint violation performance to EP-SafeBO but is more computationally efficient.

The results show that the algorithms achieve similar convergence speeds. However, the key differentiating factor is constraint violation, where SI-SafeBO outperforms both algorithms. It shows fewer violations across the 50 optimizations while maintaining consistent performance and lower variance, as indicated by the narrower box on the left of Figure 9d. Additionally, SI-SafeBO demonstrates higher computational efficiency than its counterparts, requiring less time to propose a combination of input parameters at each iteration. Further testing of the algorithms in the experimental setup of the ISOL facility is necessary.

Evaluation of Method Sensitivity

In the previous section, we showed that the proposed safe Bayesian optimization with simulation-informed Gaussian process (SI-SafeBO) for the constraints outperformed other variants of SafeBO that directly learn the classification constraints, particularly in reducing the number of constraint violations. However, in the previous experiments, it was assumed that the simulation and ground truth models coincided. Consequently, the piece-wise and true constraint functions shared the same feasibility limits. This assumption may not hold in applications where the simulation model cannot fully represent the ground truth model’s behavior. Therefore, in this section, we tested the robustness of the proposed SI-SafeBO approach through two experiments. First, we modified the simulation model relative to the ground truth model of the ISOL; second, we assumed that the simulation and ground truth models coincide, but we added noise to the design parameters of the ground truth model (quadrupole voltages $\mathbf{d}$), simulating potential hardware errors in such facilities.

Section 1 mentioned that the phase ellipse represents the radioactive ion beam (RIB), which depends on the Twiss parameters ($\alpha ,\beta ,\gamma$). A change in the initial Twiss parameters alters the extracted RIB, resulting in different beam envelopes and sizes for the same quadrupole voltages $\mathbf{d}$ compared to a RIB that maintains the original initial Twiss parameters. Consequently, the objective function and constraints differ, making the ground truth model of the ISOL and the simulation model distinct, as required for the first experiment. Therefore, we use the original Twiss parameters ${\left[\begin{array}{c}{\beta }_{x,y}\left(0\right),{\alpha }_{x,y}\left(0\right),{\gamma }_{x,y}\left(0\right)\end{array}\right]}^T={\left[\begin{array}{c}0.323,-1.233,7.809\end{array}\right]}^T$, which were used during the design of the ISOL@MYRRHA transport beamline, as the baseline. The modified Twiss parameters are generated using a normal distribution, with the baseline values as the mean and a standard deviation of $0.15$. Then, we select two sets of Twiss parameters from ten randomly generated sets. One is with the shortest Euclidean distance to the baseline parameters, and one is with the largest distance. As the Laplace approximation and expectation propagation directly learn the classification constraint, we modify the Twiss parameters of the ground truth ISOL model and fix the simulation model parameters. This setup introduces the difference between the simulation and ground truth models, enabling a comparison with the algorithms that directly learn classification constraints.

The set of Twiss parameters closest to the baseline is denoted as ${\chi }_{\mathrm{b}}$, and ${\chi }_{\mathrm{w}}$ denotes the set with the largest distance. Table 2 presents their numerical values, and Figure 10 depicts the phase ellipse of the extracted RIBs. We assume the beam has the same initial Twiss parameters in both directions, so the phase ellipse is the same in the horizontal and vertical directions. We ran the algorithm 50 times for this experiment. For the Twiss parameters ${\chi }_{\mathrm{b}}$, the algorithm started with the same 10 initial quadrupole voltages from Section 3. In contrast, for the Twiss parameters ${\chi }_{\mathrm{w}}$, certain initial locations were adjusted. This modification was required because the objective function value for some initial locations with these Twiss parameters was too close to the objective threshold $f^{\mathrm{trh}}$. We compared SI-SafeBO with LA-SafeBO and EP-SafeBO. Figure 11a,b and Figure 11c,d present the results for convergence for the Twiss parameters ${\chi }_{\mathrm{b}}$ and ${\chi }_{\mathrm{w}}$, respectively. Discrepancies between the simulation and ground truth models do not drastically affect the convergence. SI-SafeBO achieves a convergence speed similar to that of LA-SafeBO and EP-SafeBO, which directly learns the modified objective function. However, the median for LA-SafeBO and EP-SafeBO reached the objective threshold with slightly fewer objective function evaluations for the Twiss parameters ${\chi }_{\mathrm{w}}$. Nevertheless, regarding constraint violations, as depicted in Figure 12, even for ${\chi }_{\mathrm{w}}$, SI-SafeBO achieves a lower median. This suggests that even in the case of discrepancies between the simulation and ground truth models, the SI-SafeBO algorithm can find an optimal solution with fewer constraint violations than by directly learning the true classification constraints.

Moreover, there is an increasing trend in constraint violations as discrepancies between the simulation and ground truth models grow. Therefore, it is essential to note that there is a limit beyond which directly learning the constraints becomes more appropriate than using the simulation model. This is evident because the number of constraint violations for LA-SafeBO and EP-SafeBO reaches a similar median, regardless of the Twiss parameters.

Based on the results from Figure 12, it is recommended that the ISOL@MYRRHA transport beamline use the proposed SI-SafeBO as long as the Twiss parameters of the simulation model are within a normal distribution that has the ground truth parameters as the mean and a standard deviation of $0.15$. To illustrate a case where this condition is not met, we conducted an experiment using a set of Twiss parameters generated from a normal distribution with baseline values as the mean and a standard deviation of 0.25. For this experiment, we selected the Twiss parameters with the largest Euclidean distance, ${\chi }_{\mathrm{f}}={\left[\begin{array}{c}1.038,-0.952,1.835\end{array}\right]}^T$. The algorithms were run 50 times, starting from 10 different initial voltage combinations as in previous experiments. In terms of constraint violations, SI-SafeBO achieved a median of 10, while LA-SafeBO and EP-SafeBO reached medians of 9.5 and 12, respectively. Therefore, LA-SafeBO slightly outperforms SI-SafeBO. However, the 75th percentile for SI-SafeBO is 21, compared to 15 for LA-SafeBO and 17 for EP-SafeBO. Consequently, in this illustrative case, directly learning the classification constraints is advisable, as the discrepancy between the simulation model and the ground truth leads to lower performance for the SI-SafeBO algorithm compared to its counterparts.

In practice, the ground truth Twiss parameters can be estimated using diagnostic devices, such as an Allison-type scanner [51], which generates a phase space representation of the beam to facilitate Twiss parameter estimation. These devices are typically positioned after the extraction electrode, and it is common to perform a scan at the start of the tuning process to identify the extracted beam. For the simulation model’s Twiss parameters, codes like COSY INFINITY [52] allow for simulating the extracted beam and extracting its Twiss parameters. Thus, discrepancies between the ground truth and the simulation can be experimentally verified to determine the applicability of the proposed method.

For the second experiment, zero-mean Gaussian-distributed noise ${\nu }_{\mathbf{d}}$ with variance ${\sigma }_{\mathbf{d}}^2$ is added to the quadrupole voltages so that the input parameters become $\mathbf{d}\pm {\nu }_{\mathbf{d}}$. This scenario emulates hardware errors in the ISOL facility, such as offsets in the quadrupole’s power supply or errors induced by the internal controller of the power supply that does not precisely achieve the setpoint voltage. We assume that the simulation and ground truth models match, allowing us to observe only the effect of noise on the design parameters. In this experiment, we tested three different levels of noise, which are consistent with reasonable hardware errors. For instance, the largest standard deviation of the noise is ${\sigma }_{\mathbf{d}}=0.005$. It is worth noting that the quadrupole voltages are normalized, as is the noise.

The equivalent error for this noise level is approximately $\pm 35\mathrm{V}$. By introducing this level of noise to the quadrupole voltages, the Twiss parameters of the resultant beam at the entrance focal point of the dipole magnet differ from the noise-free case by approximately $10\%$, resulting in a similar proportional alteration of the final beam shape. In this experiment, we ran the algorithm 50 times, starting from the same ten different initial quadrupole voltages as in the previous experiment. Figure 13a,b depict the convergence speed and constraint violations. The convergence speed of the SI-SafeBO algorithm is similar across the three different noise levels, indicating that the algorithm is robust to changes induced by noise. Additionally, the median number of constraint violations remains the same with the introduction of noise, as shown in Figure 13d. Furthermore, the 25th and 75th percentiles are the same across the three noise levels compared to the noiseless case. These results suggest that the algorithm can effectively respond to noise in the design parameters, finding an optimal solution without significantly increasing the number of constraint violations compared to the noiseless scenario.

The results in this section show that the SI-SafeBO algorithm is a valuable approach for constrained optimization problems with classification constraints, provided a simulation model of the constraints is available. A natural progression of this work would be to further enhance the simulation-informed Gaussian process (SI-GP). One possible improvement would be to increase the robustness and reliability of the SI-GP by efficiently leveraging observations from different sources through multi-fidelity GP [53]. This approach would enable the algorithm to incorporate data from computationally inexpensive low-fidelity simulations alongside costly observations of the ground truth system. Another important direction to better utilize available information is to draw on insights from the emerging field of emulation by using a GP emulator [54] for the simulation model. This approach is especially useful in applications where the simulation model is computationally expensive to run and would allow for the identification of potentially safe regions in the design space prior to online optimization. One such software library is Emukit (Version 0.4.11) [55], a Python library compatible with Gaussian process frameworks like GPy (Version 1.0.7) [56].

In this paper, we tested two Gaussian process classifiers for directly learning the classification constraints: the Laplace approximation and expectation propagation. For the specific case of the ISOL@MYRRHA transport beamline, the performance of these classifiers was similar. However, this may not hold in other applications, where different alternatives could be advantageous. For instance, variational GP is well suited for capturing complex posterior distributions by minimizing discrepancies between the approximated and exact posterior distributions [24]. Additionally, when selecting a specific GP classifier for online optimization, it is important to consider the trade-off between classification accuracy and computational efficiency.

4. Conclusions

This paper introduces safe Bayesian optimization with a simulation-informed Gaussian process for the constraints (SI-SafeBO), a method designed to facilitate safe Bayesian optimization in scenarios involving classification constraints. SI-SafeBO requires only the availability of a simulation model of the constraints. It replaces the true classification constraint with a piece-wise function that can be learned using a Gaussian process (GP), thus enabling the application of safe Bayesian optimization. The piece-wise constraint is constructed by combining observations from the simulation model’s constraint with those from the true classification constraint. To validate SI-SafeBO’s effectiveness, we applied it to a computational model of the ISOL@MYRRHA RIB line. This study compares SI-SafeBO with safe Bayesian optimization using Laplace approximation for the constraints (LA-SafeBO) and safe Bayesian optimization with expectation propagation (EP-SafeBO). LA-SafeBO and EP-SafeBO are based on a modification of safe Bayesian optimization, enabling SafeBO to effectively use a GP classifier that learns the true classification constraints through Laplace approximation and expectation propagation.

The results demonstrate that SI-SafeBO outperformed LA-SafeBO and EP-SafeBO by reducing constraint violations by at least $80\%$ compared to these methods that directly learn the classification constraint. SI-SafeBO effectively learns the feasible design space limits with minimal constraint violations. Furthermore, sensitivity analyses were performed regarding discrepancies between the simulation model and the ground truth system and noise in the design parameters. The results showed that the method copes well with discrepancies between the simulation model and the ground truth system, as long as the Twiss parameters of the simulation model fall within a normal distribution centered on the Twiss parameters of the ground truth system, with a standard deviation of 0.15. Beyond these limits, learning the true classification constraints is advisable. With its significant reduction in constraint violations, SI-SafeBO proves to be well suited for applications where optimal solutions are required without frequent safety interlock triggers.

Author Contributions

Conceptualization, S.R.G., I.D.B. and S.D.; methodology, S.R.G.; software, S.R.G.; validation, S.R.G.; formal analysis, S.R.G., I.D.B., J.P.R. and M.D.; investigation, S.R.G.; resources, J.P.R., M.D. and L.P.; data curation, S.R.G.; writing—original draft preparation, S.R.G.; writing—review and editing, S.D., J.P.R., M.D., I.D.B. and L.P.; visualization, S.R.G.; supervision, S.D., M.D., J.P.R. and L.P.; project administration, L.P.; funding acquisition, L.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors express their gratitude to the SCK CEN Academy for funding a PhD fellowship. The work was also co-funded by the MYRRHA project at SCK CEN.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

BO	Bayesian optimization
CERN	European Organization for Nuclear Research
EP-SafeBO	Safe Bayesian optimization with expectation propagation for the constraints
GP	Gaussian process
ISAC	Isotope separator and accelerator
ISOL	Isotope separation online
ISOLDE	Isotope separation online device
LA-SafeBO	Safe Bayesian optimization with Laplace approximation for the constraints
LINAC	Linear accelerator
MRP	Mass resolving power
MYRRHA	Multipurpose Hybrid Research Reactor for High-Tech Applications
SafeBO	Safe Bayesian optimization
PSVM	Probabilistic support vector machines
RIB	Radioactive ion beam
SCK CEN	Belgian Nuclear Research Centre
SI-GP	Simulation-informed Gaussian process
SI-SafeBO	Safe Bayesian optimization with simulation-informed Gaussian process for the constraints
SVM	Support vector machines
TISA	Target and ion source assembly
TRIUMF	TRI-University Meson Facility

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic of the safe Bayesian optimization with classification constraints. (a) An approach that directly learns the true classification constraint. (b) An approach that uses the simulation-informed GP for the constraints. GP refers to the Gaussian process, SI to simulation-informed, LA to Laplace approximation, and EP to expectation propagation. The design parameters selected by SafeBO at each iteration are denoted as [Forumla omitted. See PDF.].

Figure 2. Graphical representation of an isotope separation online (ISOL) system. The shaded area represents the extracted radioactive ion beam (RIB) particles, and the thick line represents their [Forumla omitted. See PDF.] envelope. Only the RIB in the horizontal direction (plane [Forumla omitted. See PDF.]) is shown for ease of visualization.

Figure 3. Graphical representation of the phase ellipse. [Forumla omitted. See PDF.] is a general representation of the horizontal x or vertical direction y.

Figure 4. Safe Bayesian optimization illustrated with a one-dimensional variable d from the set [Forumla omitted. See PDF.]. (a) Construction of the safe set [Forumla omitted. See PDF.] using the GP’s mean [Forumla omitted. See PDF.] and standard deviation [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] as the actual constraint. [Forumla omitted. See PDF.] denotes constraint observations, and [Forumla omitted. See PDF.] is the search subspace. (b) Construction of the minimizers’ set [Forumla omitted. See PDF.] using the GP’s mean [Forumla omitted. See PDF.] and standard deviation [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] as the actual objective. [Forumla omitted. See PDF.] are the GP’s lower and upper uncertainty bounds. [Forumla omitted. See PDF.] represents objective function observations. (c) Construction of the expanders’ set [Forumla omitted. See PDF.], with the artificial location denoted as [Forumla omitted. See PDF.] and the updated GP’s mean and variance as with the artificial location ([Forumla omitted. See PDF.], [Forumla omitted. See PDF.]) (d) Selection of the next best location [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.] representing the uncertainty bounds width.

Figure 5. Graphical representation of the logistic [Forumla omitted. See PDF.] and probit functions [Forumla omitted. See PDF.], where [Forumla omitted. See PDF.] is the cumulative distribution function of the normal distribution.

Figure 6. An illustrative example of constraint modeling using the Laplace approximation. (a) Before a constraint violation occurs, the probabilistic classification ([Forumla omitted. See PDF.]) identifies the entire design space as feasible. However, the suggested location [Forumla omitted. See PDF.] lies in the unsafe region because the safe set [Forumla omitted. See PDF.] includes these locations, as the predicted upper bound [Forumla omitted. See PDF.] falls below the safety threshold (−0.5) even in unsafe areas. (b) After 30 iterations, the probabilistic classification [Forumla omitted. See PDF.] closely aligns with the true classification constraint function [Forumla omitted. See PDF.], though this accuracy is achieved only after violating the constraints four times. The set of observations is denoted as [Forumla omitted. See PDF.]. The terms [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] refer to the search subspace, the set of minimizers, and the set of expanders.

Figure 7. An illustrative example of constraint modeling using the simulation-informed Gaussian process: (a) after two iterations, (b) after 13 iterations, (c) after 30 iterations, and (d) when constraint violations occur. [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] represent the GP’s mean and standard deviation. [Forumla omitted. See PDF.] refers to the set of observations. [Forumla omitted. See PDF.] denotes the proposed piece-wise constraint function, [Forumla omitted. See PDF.] denotes the suggested location for the next iteration, and [Forumla omitted. See PDF.] is the safety threshold.

Figure 8. An illustrative example of constraint modeling using the simulation-informed Gaussian process (SI-GP) when the simulation model [Forumla omitted. See PDF.] differs from the true classification constraint [Forumla omitted. See PDF.]. (a) The selection engine considers the mean of the SI-GP at the suggested location [Forumla omitted. See PDF.] when the constraint is infeasible in the simulation model [Forumla omitted. See PDF.] but feasible in the true constraint. (b) SI-GP constraint model after 30 iterations, showing correct learning of the feasible design space. [Forumla omitted. See PDF.] refers to the set of observations.

Figure 9. Results of ISOL beamline optimization. (a) Convergence of the SI-SafeBO, (b) LA-SafeBO, and (c) EP-SafeBO algorithms. Thick lines represent the median of the minimum objective function value f from Equation (1) across 50 independent optimizations. (d) Variation in constraint violations across the 50 optimizations, where the color coding corresponds to the three different algorithms presented in panels (a–c).

Figure 10. Phase ellipses for varying Twiss parameters.

Figure 11. Convergence plot of ISOL transport beamline optimization for varying Twiss parameters. (a,b) [Forumla omitted. See PDF.] corresponding to the shortest and (c,d) [Forumla omitted. See PDF.] to the largest Euclidean distance from the baseline Twiss parameters. Thick lines indicate the median of the minimum objective function value f from Equation (1) averaged over 50 independent optimization runs. SI-SafeBO is repeated in the right plots for ease of comparison.

Figure 12. Variation of the constraint violation for ISOL transport beamline optimization with different Twiss parameters. [Forumla omitted. See PDF.] corresponds to the shortest and [Forumla omitted. See PDF.] to the largest Euclidean distance from the baseline Twiss parameters. The red plus sign refers to an outlier.

Figure 13. Results of ISOL beamline optimization for different standard deviations of input parameter noise [Forumla omitted. See PDF.]. (a) Convergence of the SI-SafeBO algorithm for (a) [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.]. Thick lines indicate the median of the minimum objective function value f from Equation (1) across 50 independent optimizations. (e) Variation in the constraint violations across the 50 optimizations, where the color coding corresponds to the different standard deviations of input parameter noise, as presented in panels (a–d).

Table 1

Comparison of time to execute tested algorithms.

Algorithm	Time to Execute SafeBO (s)
	Mean	Max
SI-SafeBO	0.61	1.27
LA-SafeBO	0.70	1.69
EP-SafeBO	1.15	3.17

Table 2

Twiss parameters of the ISOL transport beamline.

Twiss Parameters $\left[\begin{array}{c}{\mathit{\beta }}_{\mathit{x},\mathit{y}},{\mathit{\alpha }}_{\mathit{x},\mathit{y}},{\mathit{\gamma }}_{\mathit{x},\mathit{y}}\end{array}\right]$
Baseline	${\mathbf{\chi }}_{\mathbf{b}}$	${\mathbf{\chi }}_{\mathbf{w}}$
${\left[\begin{array}{c}0.323,-1.233,7.809\end{array}\right]}^T$	${\left[\begin{array}{c}0.329,-1.182,7.279\end{array}\right]}^T$	${\left[\begin{array}{c}0.633,-1.665,5.953\end{array}\right]}^T$

Appendix A

This section introduces the procedure for selecting the parameters ${\beta }_x^{\mathrm{des}}(z=B),{\alpha }_y^{\mathrm{des}}(z=B)$, and${\beta }_y^{\mathrm{des}}(z=B)$ from Equation (1). The dipole magnet at ISOL@MYRRHA has a pole face angle $\psi =31.5^{\circ }$, a bending radius $\rho =0.416\mathrm{m}$, a length ${\mathrm{L}}_m=0.7241\mathrm{m}$, a dispersion coefficient $D_m=1.191$, and a bending angle $\theta ={90}^{\circ }$. Section 2.1 mentioned that the beam must reach a specific size in the horizontal direction at the entrance focal point of the dipole magnet to achieve a specific mass resolving power (MRP) and thereby separate contaminants up to a specific mass resolution. At ISOL@MYRRHA, a mass resolving power (MRP) of 300 serves as the baseline for the design of the transport beamline. Additionally, an emittance value of $ϵ=20\pi$ mm mrad is selected for the baseline design. The MRP and the beam size relate through the expression $\mathrm{MRP}=D_m/\left[2{\sigma }_x^{\mathrm{des}}(z=B)\right]$ [47] so that the required MRP is met with a ${\beta }_x^{\mathrm{des}}(z=B)={\left[{\sigma }_x^{\mathrm{des}}(z=B)\right]}^2/ϵ\le 0.2$. It is worth noting that the higher the MRP, the better, as a higher resolution can be achieved, and thus, a purer beam can be extracted. Therefore, ${\sigma }_x^{\mathrm{des}}(z=B)=0$ is the ideal choice. However, in this case, we want to attain an MRP near 300, and therefore, we selected ${\beta }_x^{\mathrm{des}}(z=B)=0.13$.

ISOL’s transport beamline is modeled through a matrix transformation that computes the Twiss parameters at any location z along the reference trajectory. Given an initial set of parameters at $z=0$, the transformation is defined as (A1)$$\left[\begin{array}{c}\beta \left(z\right)\\ \alpha \left(z\right)\\ \gamma \left(z\right)\end{array}\right]=\left[\begin{array}{ccc}{C}^2& -2CS& S^2\\ -C{C}^{\prime }& C{S}^{\prime }+C^{\prime }S& -S{S}^{\prime }\\ C^{\prime 2}& -2C^{\prime }{S}^{\prime }& S^{\prime 2}\end{array}\right]\left[\begin{array}{c}\beta \left(0\right)\\ \alpha \left(0\right)\\ \gamma \left(0\right)\end{array}\right],$$ where the notation inside the brackets directly specifies the location z at which the Twiss parameters are obtained. The terms $C,S,C^{\prime },S^{\prime }$ inside the matrix from Equation (A1) are computed as in Equation (A2).

Satisfying the condition ${\left[\begin{array}{c}{\beta }_y^{\mathrm{des}}\left(J\right),{\alpha }_y^{\mathrm{des}}\left(J\right),{\gamma }_y^{\mathrm{des}}\left(J\right)\end{array}\right]}^T={\left[\begin{array}{c}{\beta }_y^{\mathrm{des}}\left(B\right),-{\alpha }_y^{\mathrm{des}}\left(B\right),{\gamma }_y^{\mathrm{des}}\left(B\right)\end{array}\right]}^T$ results in opposite angles at the entry and exit focal points of the dipole magnet in the vertical direction. This condition can be met by applying it to the matrix transformation from Equation (A1) in the vertical direction from point B to point J in Figure 2a. Specifically, ${\left[\begin{array}{c}\beta \left(0\right),\alpha \left(0\right),\gamma \left(0\right)\end{array}\right]}^T$ is replaced with ${\left[\begin{array}{c}{\beta }_y^{\mathrm{des}}\left(B\right),-{\alpha }_y^{\mathrm{des}}\left(B\right),{\gamma }_y^{\mathrm{des}}\left(B\right)\end{array}\right]}^T$, resulting in the matrix presented in Equation (A3). The relation $[{\gamma }_y^{\mathrm{des}}\left(B\right)=1+{\left[{\alpha }_y^{\mathrm{des}}\left(B\right)\right]}^2]/{\beta }_y^{\mathrm{des}}\left(B\right)$ [46] is adopted to eliminate ${\gamma }_y^{\mathrm{des}}\left(B\right)$ from Equation (A3), leading to the system of equations in (A4). Equation (A5) provides the solution to this system of equations for the dipole magnet at ISOL@MYRRHA. (A2)$$\begin{array}{ccc}\hfill C& =& {\displaystyle \frac{{\rho }^2-2\rho Stan\psi +{\mathrm{L}}_{m}S{(tan\psi )}^2-{\mathrm{L}}_m\rho tan\psi }{{\rho }^2}},\hfill \\ \hfill S& =& {\displaystyle \frac{(\rho -Stan\psi )({\mathrm{L}}_m\rho +2\rho S-{\mathrm{L}}_{m}Stan\psi )}{{\rho }^2}},\hfill \\ \hfill C^{\prime }& =& -{\displaystyle \frac{tan\psi (2\rho -{\mathrm{L}}_{m}tan\psi )}{{\rho }^2}},\hfill \\ \hfill S^{\prime }& =& {\displaystyle \frac{{\rho }^2-2\rho Stan\psi +{\mathrm{L}}_{m}S{(tan\psi )}^2-{\mathrm{L}}_m\rho tan\psi }{{\rho }^2}}.\hfill \end{array}$$ (A3)$$\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]=\left[\begin{array}{ccc}{C}^2-1& -2CS& S^2\\ -C{C}^{\prime }& C{S}^{\prime }+C^{\prime }S+1& -S{S}^{\prime }\\ C^{\prime 2}& -2C^{\prime }{S}^{\prime }& S^{\prime 2}-1\end{array}\right]\left[\begin{array}{c}{\beta }_y^{\mathrm{des}}\left(B\right)\\ {\alpha }_y^{\mathrm{des}}\left(B\right)\\ {\gamma }_y^{\mathrm{des}}\left(B\right)\end{array}\right].$$ (A4)$$\left[\begin{array}{c}-S^2\\ S{S}^{\prime }\end{array}\right]=\left[\begin{array}{ccc}{C}^2-1& -2CS& S^2\\ -C{C}^{\prime }& C{S}^{\prime }+C^{\prime }S+1& -S{S}^{\prime }\end{array}\right]\times \left[\begin{array}{c}{\left[{\beta }_y^{\mathrm{des}}\left(B\right)\right]}^2\\ {\alpha }_y^{\mathrm{des}}\left(B\right){\beta }_y^{\mathrm{des}}\left(B\right)\\ {\left[{\alpha }_y^{\mathrm{des}}\left(B\right)\right]}^2\end{array}\right].$$ (A5)$$\begin{array}{c}\hfill \left[\begin{array}{c}{\alpha }_y^{\mathrm{des}}\left(B\right)=c_4{\beta }_y^{\mathrm{des}}\left(B\right)\\ {\left[{\alpha }_y^{\mathrm{des}}\left(B\right)\right]}^2=c_5{\left[{\beta }_y^{\mathrm{des}}\left(B\right)\right]}^2-1\end{array}\right]\\ \hfill {\beta }_y^{\mathrm{des}}\left(B\right)=\sqrt{{\displaystyle \frac{1}{(c_5-c_4^2)}}}=1.16\mathrm{mm},\\ \hfill {\alpha }_y^{\mathrm{des}}\left(B\right)=c_4\sqrt{{\displaystyle \frac{1}{(c_5-c_4^2)}}}=1.43\mathrm{rad}.\end{array}$$

Finally, the weighting factors $c_1,c_2,\mathrm{and}{c}_3$ from Equation (1) were selected based on trial and error with multiple simulations to achieve a good trade-off between obtaining the desired MRP and the periodic solution in the vertical direction at the dipole magnet. However, higher priority was given to the beam size, resulting in a higher weight, $c_1=0.7$, compared to the weights for the periodic solution $c_2,c_3=0.15$.

Appendix B

This section contains the results of three experiments performed to select the appropriate parameters ($\kappa ,\tau ,\eta$) crucial for the SI-SafeBO algorithm. The algorithm was run ten times for each parameter tested, starting from the same initial locations from Section 2.3.2. The selection criteria are the number of constraint violations and the number of objective function evaluations from Equation (1). For the latter, the median of the ten runs is computed, and the convergence is designated as the number of objective function evaluations that the median takes to reach the objective function threshold $f^{\mathrm{trh}}.$ Table A1 presents the results for $\kappa$, where the value $5.5$ best achieves the lowest constraint violation and convergence. Regarding $\tau$, its value is expressed as a percentage of the true safety threshold. Interestingly, for the case of the ISOL@MYRRHA, not adding the safety factor achieves the lowest constraint violations and the quickest convergence. As the value of $\tau$ increases, the algorithm tends to take longer to converge due to the conservatism in exploring the design space. Finally, the selected value for $\eta$ was $0.1$ because it achieves a better trade-off between constraint violations and convergence than $0.05$, especially in convergence, and the constraint violation difference is slight.

Table A1

Constraint violations and convergence dependency on $\kappa$ for 10 runs of SI-SafeBO, where the highlighted value represents the selected one for the experiments in Section 3.

$\mathit{\kappa }$	Constraint Violations			Convergence(Obj. Evaluations)
	Median	25th Percentile	75th Percentile
2	8	6	11	33
2.5	8.5	7	10	48
3	5.5	5	9	39
3.5	4.5	3	5	51
4	3	2	5	36
4.5	3	2	4	27
5	2	1	5	48
5.5	2	1	3	26
6	2.5	1	4	48

Table A2

Constraint violations and convergence dependency on $\tau$ for 10 runs of SI-SafeBO, where the highlighted value represents the selected one for the experiments in Section 3.

$\mathit{\tau }(\%)$	Constraint Violations			Convergence(Obj. Evaluations)
	Median	25th Percentile	75th Percentile
0	2	1	3	26
2.5	3	1	4	41
5	2	2	4	36
7.5	2.5	1	4	47
10	3	1	5	44

Table A3

Constraint violations and convergence dependency on $\eta$ for 10 runs of SI-SafeBO, where the highlighted value represents the selected one for the experiments in Section 3.

$\mathit{\eta }$	Constraint Violations			Convergence(Obj. Evaluations)
	Median	25th Percentile	75th Percentile
0.05	1	1	3	31
0.1	2	1	3	26
0.15	3	2	4	41
0.2	2	2	4	39
0.25	3	2	5	37

References

1. Gui, Y.; Zhan, D.; Li, T. Taking another step: A simple approach to high-dimensional Bayesian optimization. Inf. Sci.; 2024; 679, 121056. [DOI: https://dx.doi.org/10.1016/j.ins.2024.121056]

2. Jablonka, K.M.; Jothiappan, G.M.; Wang, S.; Smit, B.; Yoo, B. Bias free multiobjective active learning for materials design and discovery. Nat. Commun.; 2021; 12, 2312. [DOI: https://dx.doi.org/10.1038/s41467-021-22437-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33875649]

3. Hanuka, A.; Huang, X.; Shtalenkova, J.; Kennedy, D.; Edelen, A.; Zhang, Z.; Lalchand, V.R.; Ratner, D.; Duris, J. Physics model-informed Gaussian process for online optimization of particle accelerators. Phys. Rev. Accel. Beams; 2021; 24, 72802. [DOI: https://dx.doi.org/10.1103/PhysRevAccelBeams.24.072802]

4. Ferran Pousa, A.; Jalas, S.; Kirchen, M.; Martinez de la Ossa, A.; Thévenet, M.; Hudson, S.; Larson, J.; Huebl, A.; Vay, J.L.; Lehe, R. Bayesian optimization of laser-plasma accelerators assisted by reduced physical models. Phys. Rev. Accel. Beams; 2023; 26, 084601. [DOI: https://dx.doi.org/10.1103/PhysRevAccelBeams.26.084601]

5. Zhang, Z.; Song, M.; Huang, X. Optimization method to compensate accelerator performance drifts. Phys. Rev. Accel. Beams; 2022; 25, 122801. [DOI: https://dx.doi.org/10.1103/PhysRevAccelBeams.25.122801]

6. Awal, A.; Hetzel, J.; Gebel, R.; Kamerdzhiev, V.; Pretz, J. Optimization of the injection beam line at the Cooler Synchrotron COSY using Bayesian Optimization. J. Instrum.; 2023; 18, P04010. [DOI: https://dx.doi.org/10.1088/1748-0221/18/04/P04010]

7. Morita, Y.; Washio, T.; Nakashima, Y. Accelerator tuning method using autoencoder and Bayesian optimization. Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrometers Detect. Assoc. Equip.; 2023; 1057, 168730. [DOI: https://dx.doi.org/10.1016/j.nima.2023.168730]

8. Kaiser, J.; Xu, C.; Eichler, A.; Santamaria Garcia, A.; Stein, O.; Bründermann, E.; Kuropka, W.; Dinter, H.; Mayet, F.; Vinatier, T. et al. Reinforcement learning-trained optimisers and Bayesian optimisation for online particle accelerator tuning. Sci. Rep.; 2024; 14, 15733. [DOI: https://dx.doi.org/10.1038/s41598-024-66263-y]

9. Shahriari, B.; Swersky, K.; Wang, Z.; Adams, R.P.; de Freitas, N. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proc. IEEE; 2016; 104, pp. 148-175. [DOI: https://dx.doi.org/10.1109/JPROC.2015.2494218]

10. Jones, D.R.; Schonlau, M.; Welch, W.J. Efficient Global Optimization of Expensive Black-Box Functions. J. Glob. Optim.; 1998; 13, pp. 455-492. [DOI: https://dx.doi.org/10.1023/A:1008306431147]

11. Vinter-Hviid, F.; Sloth, C.; Savarimuthu, T.R.; Iturrate, I. Safe contact-based robot active search using Bayesian optimization and control barrier functions. Front. Robot. AI; 2024; 11, 1344367. [DOI: https://dx.doi.org/10.3389/frobt.2024.1344367] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38741717]

12. Berkenkamp, F.; Krause, A.; Schoellig, A.P. Bayesian optimization with safety constraints: Safe and automatic parameter tuning in robotics. Mach. Learn.; 2023; 112, pp. 3713-3747. [DOI: https://dx.doi.org/10.1007/s10994-021-06019-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37692295]

13. Cole, E.R.; Connolly, M.J.; Ghetiya, M.; Sendi, M.E.S.; Kashlan, A.; Eggers, T.E.; Gross, R.E. SAFE-OPT: A Bayesian optimization algorithm for learning optimal deep brain stimulation parameters with safety constraints. bioRxiv; 2024; [DOI: https://dx.doi.org/10.1088/1741-2552/ad6cf3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39116891]

14. Kirschner, J.; Mutný, M.; Krause, A.; Coello de Portugal, J.; Hiller, N.; Snuverink, J. Tuning particle accelerators with safety constraints using Bayesian optimization. Phys. Rev. Accel. Beams; 2022; 25, 062802. [DOI: https://dx.doi.org/10.1103/PhysRevAccelBeams.25.062802]

15. Gelbart, M.A.; Snoek, J.; Adams, R.P. Bayesian Optimization with Unknown Constraints. arXiv; 2014; [DOI: https://dx.doi.org/10.48550/arXiv.1403.5607] arXiv: 1403.5607

16. Gardner, J.R.; Kusner, M.J.; Xu, Z.E.; Weinberger, K.Q.; Cunningham, J.P. Bayesian optimization with inequality constraints. Proceedings of the ICML; Beijing, China, 21–26 June 2014; Volume 2014, pp. 937-945.

17. Tfaily, A.; Diouane, Y.; Bartoli, N.; Kokkolaras, M. Bayesian optimization with hidden constraints for aircraft design. Struct. Multidiscip. Optim.; 2024; 67, 123. [DOI: https://dx.doi.org/10.1007/s00158-024-03833-8]

18. Perrone, V.; Shcherbatyi, I.; Jenatton, R.; Archambeau, C.; Seeger, M. Constrained Bayesian optimization with max-value entropy search. Proceedings of the NeurIPS 2019 Workshop on Metalearning; Vancouver, BC, Canada, 13–14 December 2019.

19. Cunningham, P.; Delany, S.J. k-Nearest Neighbour Classifiers—A Tutorial. ACM Comput. Surv.; 2021; 54, 128. [DOI: https://dx.doi.org/10.1145/3459665]

20. Kotsiantis, S.B. Decision trees: A recent overview. Artif. Intell. Rev.; 2013; 39, pp. 261-283. [DOI: https://dx.doi.org/10.1007/s10462-011-9272-4]

21. Platt, J. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Adv. Large Margin Classif.; 1999; 10, pp. 61-74.

22. Aggarwal, C.C. Data Mining: The Textbook; Springer: Berlin/Heidelberg, Germany, 2015; Volume 1, pp. 310-312.

23. Basudhar, A.; Dribusch, C.; Lacaze, S.; Missoum, S. Constrained efficient global optimization with support vector machines. Struct. Multidiscip. Optim.; 2012; 46, pp. 201-221. [DOI: https://dx.doi.org/10.1007/s00158-011-0745-5]

24. Nickisch, H.; Rasmussen, C.E. Approximations for Binary Gaussian Process Classification. J. Mach. Learn. Res.; 2008; 9, pp. 2035-2078.

25. Zhong, M.; Lotte, F.; Girolami, M.; Lécuyer, A. Classifying EEG for brain computer interfaces using Gaussian processes. Pattern Recognit. Lett.; 2008; 29, pp. 354-359. [DOI: https://dx.doi.org/10.1016/j.patrec.2007.10.009]

26. Zhang, N.; Xiong, J.; Zhong, J.; Leatham, K. Gaussian Process Regression Method for Classification for High-Dimensional Data with Limited Samples. Proceedings of the 2018 Eighth International Conference on Information Science and Technology (ICIST); Cordoba, Spain; Granada, Spain; Seville, Spain, 30 June–6 July 2018; pp. 358-363. [DOI: https://dx.doi.org/10.1109/ICIST.2018.8426077]

27. Sui, Y.; Gotovos, A.; Burdick, J.; Krause, A. Safe Exploration for Optimization with Gaussian Processes. Proceedings of the 32nd International Conference on Machine Learning; Lille, France, 6–11 July 2015; Volume 37, pp. 997-1005.

28. Kirschner, J.; Mutny, M.; Hiller, N.; Ischebeck, R.; Krause, A. Adaptive and Safe Bayesian Optimization in High Dimensions via One-Dimensional Subspaces. Proceedings of the 36th International Conference on Machine Learning; Long Beach, CA, USA, 9–15 June 2019; Volume 97, pp. 3429-3438.

29. Roussel, R.; Hanuka, A.; Edelen, A. Multiobjective Bayesian optimization for online accelerator tuning. Phys. Rev. Accel. Beams; 2021; 24, 62801. [DOI: https://dx.doi.org/10.1103/PhysRevAccelBeams.24.062801]

30. Berkenkamp, F.; Schoellig, A.P.; Krause, A. Safe controller optimization for quadrotors with Gaussian processes. Proceedings of the 2016 IEEE International Conference on Robotics and Automation (ICRA); Stockholm, Sweden, 16–21 May 2016; pp. 491-496.

31. Sui, Y.; Zhuang, V.; Burdick, J.; Yue, Y. Stagewise Safe Bayesian Optimization with Gaussian Processes. Proceedings of the 35th International Conference on Machine Learning; Stockholm, Sweden, 10–15 July 2018; Volume 80, pp. 4781-4789.

32. Krishnamoorthy, D.; Doyle, F.J. Safe Bayesian Optimization Using Interior-Point Methods—Applied to Personalized Insulin Dose Guidance. IEEE Control Syst. Lett.; 2022; 6, pp. 2834-2839. [DOI: https://dx.doi.org/10.1109/LCSYS.2022.3179330]

33. Chan, K.J.; Paulson, J.A.; Mesbah, A. Safe Explorative Bayesian Optimization—Towards Personalized Treatments in Plasma Medicine. Proceedings of the 2023 62nd IEEE Conference on Decision and Control (CDC); Singapore, 13–15 December 2023; pp. 4106-4111. [DOI: https://dx.doi.org/10.1109/CDC49753.2023.10384190]

34. Krishnamoorthy, D.; Doyle III, F.J. Model-free real-time optimization of process systems using safe Bayesian optimization. AIChE J.; 2023; 69, e17993. [DOI: https://dx.doi.org/10.1002/aic.17993]

35. Krishnamoorthy, D.; Doyle, F.J., III. Safe and Personalized Meal Bolus Calculator for Type-1 Diabetes Using Bayesian Optimization. IEEE Trans. Biomed. Eng.; 2023; 70, pp. 1481-1492. [DOI: https://dx.doi.org/10.1109/TBME.2022.3219370]

36. Rothfuss, J.; Koenig, C.; Rupenyan, A.; Krause, A. Meta-Learning Priors for Safe Bayesian Optimization. Proceedings of the 6th Conference on Robot Learning; Auckland, New Zealand, 14–18 December 2023; Volume 205, pp. 237-265.

37. Fiedler, C.; Menn, J.; Kreisköther, L.; Trimpe, S. On Safety in Safe Bayesian Optimization. arXiv; 2024; arXiv: 2403.12948[DOI: https://dx.doi.org/10.48550/arXiv.2403.12948]

38. Zhang, Y.; Park, S.; Simeone, O. Bayesian Optimization With Formal Safety Guarantees via Online Conformal Prediction. IEEE J. Sel. Top. Signal Process.; 2024; pp. 1-15. [DOI: https://dx.doi.org/10.1109/JSTSP.2024.3422825]

39. Catherall, R.; Andreazza, W.; Breitenfeldt, M.; Dorsival, A.; Focker, G.J.; Gharsa, T.P.; Giles, T.J.; Grenard, J.L.; Locci, F.; Martins, P. et al. The ISOLDE facility. J. Phys. G Nucl. Part. Phys.; 2017; 44, 094002. [DOI: https://dx.doi.org/10.1088/1361-6471/aa7eba]

40. Ball, G.C.; Hackman, G.; Krücken, R. The TRIUMF-ISAC facility: Two decades of discovery with rare isotope beams. Phys. Scr.; 2016; 91, 93002. [DOI: https://dx.doi.org/10.1088/0031-8949/91/9/093002]

41. Rasmussen, C.E.; Williams, C.K.I. Classification. Gaussian Processes for Machine Learning; The MIT Press: Cambridge, MA, USA, 2005; pp. 33-78.

42. Camps-Valls, G.; Martino, L.; Svendsen, D.H.; Campos-Taberner, M.; Muñoz-Marí, J.; Laparra, V.; Luengo, D.; García-Haro, F.J. Physics-aware Gaussian processes in remote sensing. Appl. Soft Comput.; 2018; 68, pp. 69-82. [DOI: https://dx.doi.org/10.1016/j.asoc.2018.03.021]

43. Popescu, L.; Abderrahim, H.A.; Ooms, M.; Hasaers, K. The Belgian Nuclear Research Centre (SCK CEN). Nucl. Phys. News; 2022; 32, pp. 4-8. [DOI: https://dx.doi.org/10.1080/10619127.2022.2062995]

44. Aït Abderrahim, H.; Baeten, P.; De Bruyn, D.; Fernandez, R. MYRRHA—A multi-purpose fast spectrum research reactor. Energy Convers. Manag.; 2012; 63, pp. 4-10. [DOI: https://dx.doi.org/10.1016/j.enconman.2012.02.025]

45. Köster, U. Intense radioactive-ion beams produced with the ISOL method. Eur. Phys. J. A; 2002; 15, pp. 255-263. [DOI: https://dx.doi.org/10.1140/epja/i2001-10264-2]

46. Wiedemann, H. Particle Beams and Phase Space. Particle Accelerator Physics; Springer International Publishing: Cham, Switzerland, 2015; pp. 213-251. [DOI: https://dx.doi.org/10.1007/978-3-319-18317-6_8]

47. Maloney, J.; Baartman, R.; Marchetto, M. New design studies for TRIUMF’s ARIEL High Resolution Separator. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. At.; 2016; 376, pp. 135-139. [DOI: https://dx.doi.org/10.1016/j.nimb.2015.11.023]

48. Rasmussen, C.E.; Williams, C.K.I. Regression. Gaussian Processes for Machine Learning; The MIT Press: Cambridge, MA, USA, 2005; pp. 7-32.

49. Rasmussen, C.E.; Williams, C.K.I. Covariance Functions. Gaussian Processes for Machine Learning; The MIT Press: Cambridge, MA, USA, 2005; pp. 79-104.

50. Srinivas, N.; Krause, A.; Kakade, S.M.; Seeger, M.W. Information-Theoretic Regret Bounds for Gaussian Process Optimization in the Bandit Setting. IEEE Trans. Inf. Theory; 2012; 58, pp. 3250-3265. [DOI: https://dx.doi.org/10.1109/TIT.2011.2182033]

51. Allison, P.W.; Sherman, J.D.; Holtkamp, D.B. An Emittance Scanner for Intense Low-Energy Ion Beams. IEEE Trans. Nucl. Sci.; 1983; 30, pp. 2204-2206. [DOI: https://dx.doi.org/10.1109/TNS.1983.4332762]

52. Makino, K.; Berz, M. COSY INFINITY Version 9. Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrometers Detect. Assoc. Equip.; 2006; 558, pp. 346-350. [DOI: https://dx.doi.org/10.1016/j.nima.2005.11.109]

53. Kandasamy, K.; Dasarathy, G.; Schneider, J.; Póczos, B. Multi-fidelity Bayesian Optimisation with Continuous Approximations. Proceedings of the 34th International Conference on Machine Learning; Sydney, Australia, 6–11 August 2017; Volume 70, pp. 1799-1808.

54. Paulson, J.A.; Shao, K.; Mesbah, A. Probabilistically Robust Bayesian Optimization for Data-Driven Design of Arbitrary Controllers with Gaussian Process Emulators. Proceedings of the 2021 60th IEEE Conference on Decision and Control (CDC); Austin, TX, USA, 14–17 December 2021; pp. 3633-3639. [DOI: https://dx.doi.org/10.1109/CDC45484.2021.9683046]

55. Paleyes, A.; Mahsereci, M.; Lawrence, N.D. Emukit: A Python toolkit for decision making under uncertainty. Proceedings of the 22nd Python in Science Conference; Austin, TX, USA, 10–16 July 2023.

56. GPy. GPy: A Gaussian Process Framework in Python. 2012; Available online: http://github.com/SheffieldML/GPy (accessed on 16 August 2024).