Bayesian approach for validation of runaway electron simulations

Bayesian inference aims to establish the conditional probability distribution, $P(\mathbf{x}|D)$, called the posterior, of the uncertain input parameters, $\mathbf{x}$, given observed experimental measurements, $D$. $P(\mathbf{x}|D)$ represents the best estimate and uncertainty of the input parameters for the investigated system. Bayesian inference applies the Bayes’ theorem:

$$P(\mathbf{x}|D)=\frac{P(D|\mathbf{x})P(\mathbf{x})}{P(D)},$$

(1)

which states that the posterior is proportional to the likelihood of $D$ given $\mathbf{x}$, $P(D|\mathbf{x})$, multiplied by the prior probability distribution for $\mathbf{x}$, $P(\mathbf{x})$. The marginal probability of the experimental measurements, $P(D)$, represents an integral over all possible data generating input values, $\int P(D|\mathbf{x}^{\prime})P(\mathbf{x}^{\prime})d\mathbf{x}^{\prime}$, which would be computationally challenging to evaluate. However, it is typically sufficient to establish the relative posterior probabilities of various values of $\mathbf{x}$. Therefore, $P(D)$ does not need to be directly evaluated, and it is sufficient to apply the proportionality:

$$P(\mathbf{x}|D)\propto P(D|\mathbf{x})P(\mathbf{x}).$$

(2)

When the likelihood function, $P(D|\mathbf{x})$, is either not available analytically or cannot be evaluated within the available computational or time resources, the standard alternative is to use approximate Bayesian computation (ABC) Marin et al. (2012); Järvenpää et al. (2018). ABC aims to establish the approximate Bayesian posterior:

$$P_{\text{ABC}}(\mathbf{x}|D)\propto P(\mathbf{x})\int H(\epsilon-\Delta(D,y))P(y|\mathbf{x})dy,$$

(3)

where $y\in\mathbb{R}^{d}$ is data generated with the simulation model with input parameters $\mathbf{x}$, and $\Delta:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{>0}$ is a discrepancy function between the simulated and measured data. $\epsilon$ represents the threshold parameter controlling the trade-off between posterior estimation accuracy and efficiency, and $H(x)$ is the Heaviside step function which takes a value 1 whenever $\epsilon$ is greater than the discrepancy. Small values of $\epsilon$ lead to more accurate approximate posteriors, but also increase the computational challenge.

One of the simplest ABC algorithm that could be applied to numerically estimate the integral (Eq. 3) is rejection sampling Marin et al. (2012); Järvenpää et al. (2018); Lintusaari et al. (2017):

1. 1.

   Draw a random sample from $\mathbf{x}^{\prime}\sim P(\mathbf{x})$.

2. 2.

   Generate $y^{\prime}\sim P(y|\mathbf{x}^{\prime})$.

3. 3.

   Accept $y^{\prime}$ if $\Delta(D,y^{\prime})\leq\epsilon$.

4. 4.

   Go back to step (i) until sufficient number of accepted samples are collected.

The accepted values represent the approximate posterior distribution. The drawback of the standard rejection sampling is the number of required function evaluations. For a typical simulation tool in magnetic confinement fusion, a function evaluation takes at least several minutes and more typically hours or days. Therefore, it is not computationally feasible to collect sufficiently many samples to get an accurate ABC posterior distribution using a rejection sampler. Furthermore, a rejection sampler with a small threshold parameter, $\epsilon$, is computationally very inefficient as a large fraction of the sampled function evaluations are rejected. While it is possible to improve the efficiency by applying approaches such as Markov Chain Monte Carlo (MCMC) Marjoram et al. (2003), the direct sampling based ABC algorithms are still expected to be computationally too costly for the type of applications targeted in this work.

The inefficiency of the direct sampling based ABC approach can be circumvented by applying Bayesian optimization to traverse the space of optimized input parameters Gutmann & Corander (2016); Järvenpää et al. (2019). This leads to a probabilistic surrogate model based ABC approach that is observed to be several orders of magnitude more efficient in terms of full function evaluations than the direct sampling based ABC algorithms. At each step of the algorithm, the ABC posterior is estimated using the surrogate model as $P_{\text{ABC}}(\mathbf{x}|D)\propto P(\mathbf{x})\mathbb{P}(\Delta_{\mathbf{x}}<\epsilon)$, where the probability is computed using the probabilistic surrogate model.

Bayesian optimization (BO) offers a powerful approach for global optimization of costly-to-evaluate, non-convex functions, without access to first- or second-order derivatives Brochu et al. (2010); Shahriari et al. (2016); Frazier (2018). The problem of finding optimal values, $\mathbf{x_{*}}$, for the uncertain input parameters, $\mathbf{x}$, can be represented as a task of finding the optimum of a non-linear function $f(\mathbf{x})$ of a compact set $\mathcal{A}$, called search space in this report. If $f(\mathbf{x})$ represents the discrepancy between the synthetic and measured diagnostic data, then the problem can be formulated as:

$$\mathbf{x_{*}}=\arg\min_{\mathbf{x}\in\mathcal{A}}f(\mathbf{x}).$$

(4)

The target for a BO algorithm is to be able to traverse the search space efficiently in terms of function evaluations and to find the globally optimum solution by applying prior belief about the optimized function and by balancing exploration and exploitation of the search space Brochu et al. (2010); Shahriari et al. (2016); Frazier (2018). In exploitation, samples are collected in regions of the search space that are known to lead to near optimal function values based on prior belief, and in exploration, samples are collected in regions that encompass a large uncertainty.

A standard BO algorithm consists of two main components Brochu et al. (2010); Shahriari et al. (2016); Frazier (2018):

1. 1.

   A probabilistic model for the objective function.

2. 2.

   An acquisition function for recommending the next sampling point.

The probabilistic model represents essentially a low evaluation cost surrogate model for the objective function, and the uncertainties retained in the probabilistic model represent the degree of confidence on the surrogate model predictions. The acquisition function applies the mean and variance of the probabilistic model to balance exploitation and exploration. The collected sample values are then used to update the probabilistic surrogate model.

The usual choice for the probabilistic model is to use Gaussian Process Regression (GPR), also known as Kriging (Stein, 1999; Rasmussen & Williams, 2006, and references therein). Kriging surrogate-based optimization was previously used for parameter optimization of plasma transport codes by Rodriguez-Fernandez et al. (2018). Other examples of GPR applications in plasma physics can be found in (Ho et al., 2019; Chilenski et al., 2017, 2015; von Nessi et al., 2012; Li et al., 2013; Romero & Svensson, 2013; von Nessi & Hole, 2013, and references therein).

GPR is a Bayesian regression technique and is very powerful for interpolating small sets of data as well as retaining information about the uncertainty of the regression (Stein, 1999; Rasmussen & Williams, 2006, and references therein). Gaussian process (GP) is a stochastic process, for which any finite collection of random values has a multivariate normal distribution. GP is specified by the mean function, $m(\mathbf{x})=\mathbb{E}[f(\mathbf{x})]$, and the covariance function, $k(\mathbf{x},\mathbf{x^{\prime}})=\mathbb{E}[(f(\mathbf{x})-m(\mathbf{x}))(f(\mathbf{x^{\prime}})-m(\mathbf{x}^{\prime}))]$ Rasmussen & Williams (2006):

$$f(\mathbf{x})\sim\mathcal{GP}(m(\mathbf{x}),k(\mathbf{x},\mathbf{x^{\prime}})).$$

(5)

The GPR represents a family of functions, for which the point to point variance is described by the covariance function. Usually, the mean function is assumed as zero, although other assumptions are possible. This means that the mean of the prior assumption on variation of the objective function value when propagating from a collected data point is zero. The covariance function or kernel describes the smoothness assumption on the possible functions $f$. GP essentially describes a normal distribution over possible functions, $f$, conditioned with observations {$(\mathbf{x}_{i},f_{i})$, $i=1,...,n$}. Assuming a collection of possibly noisy observations with a Gaussian noise variance $\sigma_{n}^{2}$, the posterior probability distribution function of $f$ at point $\mathbf{x}$ is Gaussian with posterior mean $\mu_{n}(\mathbf{x})$ and posterior variance $v_{t}(\mathbf{x})+\sigma_{n}^{2}$:

$$f(\mathbf{x})\sim\mathcal{N}(\mu_{n}(\mathbf{x}),v_{t}(\mathbf{x})+\sigma_{n}^{2}).$$

(6)

Assuming a zero mean function, $m(\mathbf{x})$, the mean and variance can be obtained as

$$\mu_{n}(\mathbf{x_{*}})=\mathbf{k}_{*}^{\intercal}\mathbf{K}_{n}^{-1}\mathbf{f_{n}},$$

(7)

$$v_{n}(\mathbf{x_{*}})=k(\mathbf{x_{*}},\mathbf{x_{*}})-\mathbf{k}_{*}^{\intercal}\mathbf{K}_{n}^{-1}\mathbf{k}_{*},$$

(8)

where $\mathbf{k_{*}}$ represents the vector of covariances between the test point, $\mathbf{x_{*}}$, and the $n$ observations, $f_{n}$ is a vector of the $n$ observations, and $\mathbf{K}$ is the covariance matrix Gutmann & Corander (2016); Rasmussen & Williams (2006). Since the function evaluations in this work are deterministic, the $\sigma_{n}$ term is constrained to a low value that does not impact the predictions. An estimate of the likelihood at $w$ is given by

$$\mathbb{P}(\Delta_{w}<\epsilon)=F\left(\frac{\epsilon-\mu_{n}(w)}{\sqrt{v_{n}(w)+\sigma_{n}^{2}}}\right),$$

(9)

where $F(x)$ is the cumulative distribution function of $\mathcal{N}(0,1)$. In this work, we will set $\epsilon$ equal to the current optimal value provided by the probabilistic surrogate model, which is also the default in BOLFI Gutmann & Corander (2016).

A key step in building a GP regression is to select the covariance function or kernel Rasmussen & Williams (2006). Usually the default choice is the radial basis function (RBF), also known as squared exponential or Gaussian kernel:

$$k_{\text{RBF}}(x_{i},x_{j})=\sigma_{f}^{2}\text{exp}\left(-\sum_{i=1}^{d}\frac{\left(x_{i,k}-x_{j,k}\right)^{2}}{2l_{k}^{2}}\right),$$

(10)

where $\mathbf{l}=[l_{1},...,l_{d}]$ is a vector of covariance lengthscales for each dimension, $d$, and $\sigma_{f}^{2}$ is the variance. In the applications in this work, the single constant $\mathbf{l}$ was observed to be often too restrictive and a rational quadratic kernel (RQ) was used instead:

$$k_{\text{RQ}}(x_{i},x_{j})=\sigma_{f}^{2}\left(1+\sum_{i=1}^{d}\frac{\left(x_{i,k}-x_{j,k}\right)^{2}}{2\alpha l_{k}^{2}}\right)^{-\alpha},$$

(11)

which is equivalent to summing many RBF kernels with varying $\mathbf{l}$. The hyperparameter $\alpha$ represents the relative weighting between large and small $\mathbf{l}$ values. Before applying the model, the hyperparameters ($\mathbf{l}$, $\alpha$, $\sigma_{f}^{2}$, $\sigma_{n}$) must be optimized first. These can be estimated by maximizing the marginal log-likelihood Rasmussen & Williams (2006). The GPR library used in this work as well as in BOLFI is the Python GP framework GPy GPy (since 2012), which encompasses the applied optimization routines.

The acquisition function applies the mean and variance of the probabilistic surrogate model to recommend the next sampling point for the objective function Brochu et al. (2010); Shahriari et al. (2016); Frazier (2018). The acquisition functions are typically constructed to recommend sampling points that either encompass optimal predicted mean for the objective function, exploitation, or large uncertainty, exploration. The sampling point is selected by optimizing the acquisition function. Several acquisition functions have been developed and can be found in the reviews in Brochu et al. (2010); Shahriari et al. (2016); Frazier (2018). The acquisition function used in sequential sampling in this work is the lower confidence bound selection criterion (LCBSC), which is also the default acquisition function in BOLFI (Brochu et al., 2010; Gutmann & Corander, 2016; Srinivas et al., 2010, and references therein). This function can be written as

$$\mathcal{A}_{n}(\mathbf{x})=\mu_{n}(\mathbf{x})-\sqrt{\eta_{n}^{2}v_{n}(\mathbf{x})},$$

(12)

$$\eta_{n}^{2}=2\text{ln}\left(\frac{n^{2d+2}\pi^{2}}{3\epsilon_{\eta}}\right),$$

(13)

where the coefficient $\eta_{n}$ is the trade-off parameter between exploration and exploitation, $\epsilon_{\eta}\in(0,1)$ is a constant chosen by the user, and $d$ is the dimensionality of the search space. Optimizing the acquisition function provides a deterministic answer for the next sampling point. An example of the application of this acquisition function is shown in Section 3.2.

Since the next sampling point is obtained deterministically for a given state of the surrogate model, the approach is naturally sequential: (1) the objective function is evaluated for the sampling location provided by the optimum of the acquisition function, (2) the GPR surrogate model is updated, and (3) the acquisition function is optimized again, using the updated GPR, to recommend the next sampling point. However, with complicated, multi-dimensional optimization tasks with computationally time consuming function evaluations, it would be more attractive to conduct several objective function evaluations in parallel to each other to reduce the overall time consumption of the optimization task, especially when suitable high-performance computing (HPC) resources are available.

To conduct parallel BO, stochastic acquisition rules can be used. Various approaches for parallel BO and batch acquisition have been developed Järvenpää et al. (2018, 2019); Thompson (1933); Kandasamy et al. (2018); Chapelle & Li (2011); Palma et al. (2019). In this work, the randmaxvar approach developed by Järvenpää et al. (2019) is used as the stochastic acquisition method. The approach is based on the maxvar acquisition rule also presented in Järvenpää et al. (2019). The maxvar acquisition method recommends a sample in a location that encompasses the maximum variance of the unnormalized ABC posterior. Basically, due to the limited information, represented by the collected samples, there is uncertainty in the GPR representation of the objective function and this uncertainty is propagated as uncertainty of the unnormalized ABC posterior. The maxvar method aims to collect samples that lead to maximum reduction of this uncertainty. In the stochastic version of this method, samples are collected from the distribution that represents the variance of the unnormalized ABC posterior. Since samples are collected stochastically, several samples can be collected without updating the GPR surrogate in between the samples, enabling parallelization. Furthermore, sampling can be done asynchronously by simply updating the GPR surrogate whenever new results are added to the dictionary of collected samples and by sampling new values whenever an idle processor becomes available.

The BO and ABC approach discussed in the previous section is applied for DREAM runaway electron simulations of current quench (CQ) in the disruption of a JET discharge #95135. This was a deuterium limiter plasma with an argon massive gas injection induced disruption, described in detail by Reux et al. (2021) and Brandström (2021). DREAM is a numerical tool for self-consistently simulating the evolution of temperature, poloidal flux, and impurity densities, along with the generation and transport of REs in tokamak disruptions (Hoppe et al., 2021). The DREAM simulations in this manuscript are similar to those presented by Brandström (2021) with the exception that only a fluid model for RE electrons is used here, as kinetic simulations were not necessary for the proof-of-principle of the Bayesian approach. For a full description of the physics model in DREAM, we refer the reader to Ref. (Hoppe et al., 2021).

The simulations are started at the peak of the total plasma current, $I_{\text{p}}$, obtained during the disruption. An instantaneous thermal quench (TQ) is assumed, after which all background plasma quantities, except the electron temperature, $T_{\text{e}}$, are evolved self-consistently. The post-disruption $T_{\text{e}}$ is instead given as an uncertain input parameter. The background plasma density is obtained from pre-disruption measurement with the high resolution Thomson scattering (HRTS) (Pasqualotto et al., 2004; Frassinetti et al., 2012). Even though the uncertainty of the electron density measurement could be taken into account in the Bayesian approach, for simplicity we neglect it here.

The argon density in the plasma is obtained from the estimated amount of injected argon, volume of the vessel, and fraction of argon that is assimilated. Argon is assumed to be uniformly distributed in the plasma. While the injected amount can be obtained from the experiment, $N_{\text{Ar}}\sim 8\times 10^{20}$ atoms Brandström (2021), the assimilated fraction is given as an uncertain input parameter, $f_{\text{Ar}}$.

Hot-tail and Dreicer RE generation mechanisms were not self-consistently included in the simulations. Instead, a RE seed profile is given as an uncertain input to generate a RE beam through the avalanche mechanism. The initial total plasma current $I_{\text{p}}$ (combination of ohmic and prescribed runaway seed current) is constrained to match the experimentally measured peak value, $I_{\text{p}}\sim 1.42$ MA, by adjusting the initial electric field in order for the avalanche multiplication factor to to be constrained. The electric field is evolved self-consistently, during the CQ and RE plateau simulations. The conductivity of the wall is controlled by a characteristic wall time, $\tau_{\text{wall}}=L_{\text{ext}}/R_{\text{wall}}$, that is provided as an uncertain input parameter. Here, $L_{\text{ext}}$ is the external inductance and $R_{\text{wall}}$ the resistance of the wall.

The full list of uncertain input parameters are $T_{\text{e}}$, $f_{\text{Ar}}$, RE seed distribution, and $\tau_{\text{wall}}$. In addition, the RE seed distribution is scaled with a multiplier such that the RE plateau current matches the experimentally measured value. This is done iteratively with a binary search algorithm. The algorithm is initialized by multiplying or dividing the multiplier by $10^{3}$ at each iteration until the predicted plateau current passes the experimentally measured value, after which binary search is applied to converge to the optimum. If the multiplier is reduced below a small value, $\sim 10^{-3}$, the search algorithm stops without finding a solution that matches the RE plateau current and the objective function value at the point of reaching that threshold is propagated through the algorithm. This is simply due to the fact that it cannot be assumed that for all input values within the search space it is possible to find a multiplier that would enable matching the RE plateau current.

In the proof-of-principle examples that follow, the objective function is chosen as the L1-norm of the discrepancy between the measured and predicted $I_{\text{p}}$ during the CQ (Fig. 1):

$$\Delta=\lVert I_{\text{p}}^{\text{measured}}-I_{\text{p}}^{\text{DREAM}}\rVert_{1}.$$

(14)

Effectively this discrepancy function calculates the area between the two curves in Figure 1. To avoid accumulating excessive integration within the RE plateau, the discrepancy is only calculated between 0 and 25 ms from the current peak.

The first proof-of-principle test is to apply the Bayesian approach to search for the post TQ $T_{\text{e}}$ that minimises the discrepancy function (Eq. 14). The other uncertain input parameters are fixed as $f_{\text{Ar}}=15$%, uniform RE seed profile, and $\tau_{\text{wall}}=5$ ms. A rational quadratic kernel is used and the lengthscale is constrained to be below 1.0. This is done to avoid the model becoming overconfident in regions that have not been sampled yet. If the lengthscale converges to a large value, it can suppress exploration prematurely and prevent the algorithm from finding the optimal solution. More generally it seems that in BO algorithms it is probably better to have a model that has capability to overfit rather than underfit the data. The BO algorithm is interpolating solutions within the search space and an overfitting model will just encourage exploration while an underfitting model might not have the generalization capability to fit the solution near the optimum. The acquisition function is the LCBSC (Eq. 12) with exploration constant $\epsilon_{\eta}=0.2$. The search space for $T_{\text{e}}$ is bounded between 1.0 and 20 eV. A uniform prior uncertainty distribution is assumed.

Within less than 10 iterations, the algorithm starts to converge to the optimum value (Fig. 2). The first 3 samples are collected randomly from the uniform prior distribution. After these are collected, the GPR surrogate model is fitted to the data (Fig. 2). The mean and variance of the GPR are applied by the acquisition function to recommend the next sampling location. The algorithm proceeds by choosing the location of the minimum value of the acquisition function. By proceeding like this, the algorithm converges near the optimum value with a narrow confidence interval around the sample number 7 to 8 (Fig. 3a). Since the prior distribution is uniform, the posterior probability distribution can be obtained from the GPR by applying Eq. 9 (Fig. 3b). Finally, the simulated current value with the temperature value providing the highest posterior probability, $T_{\text{e}}\sim 5.84$ eV, can be compared to the experimental measurements (Fig. 3c). The results indicate that by only adjusting the constant post thermal quench $T_{\text{e}}$, the model is not able to reproduce the experimentally measured CQ rate. Alternatively, these results indicate that the background plasma resistivity is changing during the CQ.

After the 1D proof-of-principle, the next step is to extend the search space to include the other uncertain parameters as well. Since the RE seed is a 1D profile, it is parameterized as the probability distribution function of the gamma distribution:

$$f_{\text{RE}}(r;\alpha,\beta)=\frac{\beta^{\alpha}r^{\alpha-1}e^{-\beta r}}{\Gamma(\alpha)},$$

(15)

where $\alpha$ and $\beta$ are free parameters and $\Gamma(\alpha)$ is the gamma function (Fig. 4). The intention is to provide a general parameterization that constrains as little as possible the possible RE seed profile shapes. The search spaces for both $\alpha$ and $\beta$ are set as uniform distributions between 0.001 and 10. 10 randomly sampled RE seed profile shapes are shown in Figure 4. It can be clearly seen that this parameterization allows flexible representations of profiles peaking at the center, middle, or edge of the plasma. The search space for the argon assimilation fraction is set as uniform between 0.001% and 100%. The search space for the characteristic wall time is set such that $\text{ln}\left(\tau_{\text{wall}}/1\,\mathrm{ms}\right)$ is sampled uniformly between 0 and 7. As a result, the $\tau_{\text{wall}}$ values range between about 1 ms and 1100 ms. This allows sampling for very large $\tau_{\text{wall}}$ parameters exceeding 1.0 s, while encouraging collection of samples at low values.

Due to the expanded volume of the search space, more samples are expected to be needed than in the 1D example. Therefore, the randmaxvar acquisition function was used with a batch size of 10 samples conducted in parallel. The first 50 samples were collected by sampling the search space randomly, after which the acquisition function was used to recommend sampling locations.

Similar to the 1D search, a rational quadratic kernel is used. The kernel parameters are restricted such that the power parameter in the rational quadratic kernel was constrained to be between $10^{-10}$ and 0.03. The lengthscale constraints are altered for every batch of 10 samples. For even round batches, the lengthscales are constrained to be positive and manually initialized as the distance of the search domain for the dimension divided by the number of collected samples. After this preconditioning step, GP optimization is conducted. There is no maximum lengthscale setup during the even rounds. During the odd round batches, the lengthscales for variations are constrained to be below 1 for the $T_{\text{e}}$ and $f_{\text{Ar}}$ dimensions, below 0.5 for the $\alpha$ and $\beta$ dimensions, and below 0.1 for the $\text{ln}\left(\tau_{\text{wall}}\right)$ dimension. The lower limits for the lengthscales were set to $10^{-3}$. The even rounds perform essentially automatic relevance determination obtaining very long lengthscales for dimensions that do not show a significant impact on the objective function value, guiding the search to prioritize optimization of input parameters that have more significant impact on the objective function value. However, this approach alone would risk the surrogate model becoming overly confident early in the search and stop exploration for input parameters deemed unimportant. To counteract this risk, the odd rounds apply restrictions of lengthscales, such that the algorithm understands to explore input parameter values, which would simply be extrapolated and interpolated over by the long-lengthscale surrogate model.

The sampling algorithm finds clear optima in the search space for the background plasma temperature and characteristic wall time (Fig. 5). It can be observed that there is a global optimum at $T_{\text{e}}$ around 5 – 7 eV and $\tau_{\text{wall}}$ less than about 2.7 ms ($\text{ln}(\tau_{\text{wall}}/\text{1 ms})<1.0$), and a local optimum at $T_{\text{e}}$ about 15 – 20 eV and $\tau_{\text{wall}}$ larger than about 50 ms ($\text{ln}(\tau_{\text{wall}}/\text{1 ms})>4.0$) (Fig. 5 red and black ellipses). The two solution branches can also be observed in the plot as a function of the argon assimilation fraction, such that the low temperature solution branch reaches the optimum at values above 20 %, while the higher temperature solution branch at lower values between 5 – 30 % (Fig. 5c). The shape of the RE seed distribution does not seem to impact the discrepancy significantly (Figs. 5d, e).

After 290 samples, the algorithm estimates the global optimum to be $T_{\text{e}}\approx 6.1$ eV, $f_{\text{Ar}}\approx 54.2$ %, $\tau_{\text{wall}}\approx 1.1$ ms, $\alpha\approx 7.9$, and $\beta\approx 4.3$, where the 95 % confidence intervals for $f_{\text{Ar}}$, $\alpha$, and $\beta$ span most of the search space (Fig. 6). The local uncertainties for the global optimum can be obtained by evaluating the local properties of the approximate posterior. The approximate posterior can be extracted from the probabilistic surrogate model by applying Eq. 9. While a global ABC posterior can be obtained through, for example, MCMC sampling, the local inverse uncertainty near the global optimum is likely of more practical interest in fusion energy research and can be obtained directly from the GPR representation of the posterior. This type of analysis was done by evaluating the one dimensional posterior distribution for each search dimension from the global optimum, which can be integrated to obtain confidence intervals (Fig. 6). It should be noted that this analysis does not take into account the secondary optimum at higher temperatures as that would require non-linear analysis of the approximate posterior, while MCMC sampling of the posterior would collect some distribution in that area also. However, since multimodality of the optimized function can be observed from the discrepancy plot already (Fig. 5), it is considered more important to obtain local uncertainty estimations near the global optimum than sample the full approximate posterior. For completeness, MCMC sampling was conducted for the approximate posterior and the resulting distributions for $T_{\text{e}}$ and $\text{ln}(\tau_{\text{wall}})$ show the global optimum, also visible in (Fig. 6), as well as the secondary local optimum at higher temperatures (Fig. 7).

Analysing the convergence of the search it can be observed that after about 80 – 110 samples, the posterior distributions for $T_{\text{e}}$ and $\text{ln}(\tau_{\text{wall}})$ have converged (Figs. 6a, b). For the other parameters, the 95 % confidence intervals remain large, indicating that the discrepancy value is not very sensitive to these input parameters, as was observed in the discrepancy plot already (Fig. 5). The step observed at 200 samples highlights the stochastic nature of the GPR surrogate model fitting. Depending on the initial conditions, the optimization algorithm might find somewhat different hyperparameters for the GPR, leading to a different shape of the approximate posterior and shift of the optimum and boundaries of the confidence interval as well. However, the large steps of the optimum for $f_{\text{Ar}}$, as well as the large steps at other sample numbers for $\alpha$, and $\beta$, are a result of the relatively flat posterior distribution shapes, while posterior shape for $T_{\text{e}}$ and $\text{ln}(\tau_{\text{wall}})$ are not changed significantly.

Finally, the predicted and measured total plasma current are compared for the global optimum extracted by the algorithm (Fig. 6f). Comparing the result to the optimum case in the 1D search space example, it can be clearly seen that, with a 5D search space, the algorithm is able to obtain a significantly better fit to the experimentally measured current (Figs. 3c, 6f). In the 1D example, the current was reducing significantly faster than experimentally measured during the early part of the CQ and the end of the CQ happened several ms later than experimentally measured, such that the average L1-norm discrepancy was minimized, while the rate of change of plasma current was poorly matched (Fig. 3c). On the other hand, in the 5D example, the initial drop is also faster than measured experimentally, but soon the rate of change of the plasma current is matched to the experimentally measured rate of change, such that the end of CQ happens near the experimentally measured end of the CQ when the L1-norm between the two currents is minimized (Fig. 6f). The fact that the current is reducing faster than experimentally measured during the early parts of the CQ indicates that the plasma resistance that on average matches the plasma current evolution probably overestimates the resistance during early parts of the CQ. To address this, the final proof-of-principle test conducted in this manuscript is to allow linear variation of $T_{\text{e}}$ during the simulation by extending the search space to seven dimension.

As a next extension of the search, a parameterized variation of the background plasma $T_{e}$ during the simulation is allowed. Linear variation of $T_{e}$ is assumed, such that the search space is extended to seven dimensions, by adding final temperature, $T_{\text{e, final}}$, and the time at which $T_{\text{e, final}}$ is reached, $t_{\text{final}}$. After $t_{\text{final}}$ is reached, $T_{\text{e}}$ is assumed to stay constant. Same search spaces are used for the initial and final $T_{\text{e}}$, and the search space for $t_{\text{final}}$ is set as uniform between 1 ms and 44 ms.

The rational quadratic kernel is used in the GPR. The power of the kernel was restricted similarly to the setup in the 5D search. Similar to the 5D search, the lengthscale restrictions were altered between even and odd round batches. Batch size was set to 50. For the odd round batches, the lenghscale constraints are similar to those in the 5D search with the same lengthscale contraint applied for the initial and final $T_{e}$. For the $t_{\text{final}}$ the minimum lenghtscale is set as 10^-3 ms and maximum as 1 ms.

After 950 samples, the local 95 % confidence intervals around the optimum point, recommended by the algorithm, are $T_{\text{e, initial}}\in[8.5,11.7]$ eV, $T_{\text{e, final}}\in[3.3,5.8]$ eV, $t_{\text{final}}\in[11,19]$ ms, $\tau_{\text{wall}}\in[1.1,2.1]$ ms, $f_{\text{Ar}}\in[48,97]$ %, $\alpha\in[0.2,9.4]$, and $\beta\in[0.5,9.7]$ (Fig. 8). The optimum point is $T_{\text{e, initial}}\approx 9.5$ eV, $T_{\text{e, final}}\approx 4.3$ eV, $t_{\text{final}}\approx 14$ ms, $\tau_{\text{wall}}\approx 1.4$ ms, $f_{\text{Ar}}\approx 76$ %, $\alpha\approx 4.3$, and $\beta\approx 2.2$. With these input parameters, the predicted current quench rate is very close to the measured values (Fig. 8f). Initially, the rate is faster than measured and also the transition to runaway plateau in the simulation occurs about 1 ms earlier than measured (Fig. 8f). However, between 3 and 12 ms, the predicted current is nearly exactly on top of the measured current.

In both the 5D and 7D searches, the algorithm found optimum $\tau_{\text{wall}}$ around 1.1 - 1.4 ms. This result was somewhat surprising when the conventional prior expectation would have suggested higher values in the range of 5 to 10 ms. Since the optimization algorithm finds the set of parameters that minimize the discrepancy, it is possible that the algorithm compensates for missing physics in the model by reducing $\tau_{\text{wall}}$ below values that would actually be realistic. As a final test, further 4D and 6D optimizations were conducted, where $\tau_{\text{wall}}$ was fixed to 5 ms.

In the 4D search, the optimum point recommended by the algorithm is $T_{\text{e}}\approx 7.0$ eV, $f_{\text{Ar}}\approx 52$ %, $\alpha\approx 4.6$, and $\beta\approx 5.5$ (Fig. 9a). The local 95 % confidence intervals are $T_{\text{e}}\in[6.4,8.0]$ eV and $f_{\text{Ar}}\in[32,91]$ %, $\alpha\in[0.3,9.6]$, and $\beta\in[0.3,9.7]$. In the 6D search, the optimum point recommended by the algorithm is $T_{\text{e, initial}}\approx 16.2$ eV, $T_{\text{e, final}}\approx 5.9$ eV, $t_{\text{final}}\approx 10$ ms, $f_{\text{Ar}}\approx 82$ %, $\alpha\approx 1.3$, and $\beta\approx 6.1$ (Fig. 9b). The local 95 % confidence intervals are $T_{\text{e, initial}}\in[13.9,19]$ eV, $T_{\text{e, final}}\in[5.1,6.7]$ eV, $t_{\text{final}}\in[8,12]$ ms, $f_{\text{Ar}}\in[48,98]$ %, and $\alpha\in[0.1,8.8]$, $\beta\in[0.5,9.7]$.

As the $f_{\text{Ar}}$, $\alpha$, and $\beta$ do not impact the discrepancy significantly, the best match obtained by the 4D search seems very similar to the best match obtained by the 1D search (Figs. 3c, 9a). When allowing linearly varying background plasma $T_{\text{e}}$, the algorithm is able to find a solution that matches the experimentally measured plasma current nearly as well as in the 7D search (Figs. 8f, 9b). However, when fixing $\tau_{\text{wall}}=5$ ms, the recommended initial $T_{\text{e}}$ is increased from 9.5 to 16.2 eV, and $t_{\text{final}}$ reduced from 14 ms to 10 ms, highlighting the non-linear dependencies between the optimal input parameters.

The computational challenge of the optimization task increases with the number of dimensions of the search space. Figure 10 illustrates the discrepancy as a function of sample numbers for the 4D, 5D, 6D, and 7D search tasks in this work. Evaluating the convergence based on the sample number after which the minimum discrepancy saturates, the 4D search converges after about 40 samples, the 5D search after about 80 samples, the 6D search after about 250 samples, and the 7D after about 300 samples. Beyond this point, increasing sample numbers will reduce the uncertainty of the posterior distribution while the minimum discrepancy is not reduced anymore.

Comparing to a grid search of eight samples for each dimension, the Bayesian optimization algorithm is very efficient (Fig. 11). Beyond four dimensions, the grid search algorithm would be calling for over 10000 samples and soon become intractable. The Bayesian approach, on the other hand, obtains samples near the minimum discrepancy after a few hundred samples even in the case of the seven dimensional search space.