Automated simulation-based design via multi-fidelity active learning and optimisation for laser direct drive implosions

TL;DR

This paper presents a machine learning framework that leverages multi-fidelity simulations to efficiently optimize laser-driven inertial fusion designs, reducing computational costs while ensuring robustness against hydrodynamic instabilities.

Contribution

It introduces an ensemble of neural network surrogates trained on 1D and 2D data for robust design optimization in inertial fusion experiments.

Findings

Surrogate models effectively predict robust designs.

Optimized designs achieve high gain in 2D simulations.

Framework reduces reliance on costly multi-dimensional simulations.

Abstract

The design of inertial fusion experiments is a complex task as driver energy must be delivered in a precise manner to a structured target to achieve a fast, but hydrodynamically stable, implosion. Radiation-hydrodynamics simulation codes are an essential tool in this design process. However, multi-dimensional simulations that capture hydrodynamic instabilities are more computationally expensive than optimistic, 1D, spherically symmetric simulations which are often the primary design tool. In this work, we develop a machine learning framework that aims to effectively use information from a large number of 1D simulations to inform design in the presence of hydrodynamic instabilities. We use an ensemble of neural network surrogate models trained on both 1D and 2D data to capture the space of good designs, i.e. those that are robust to hydrodynamic instabilities. We use this surrogate to perform Bayesian optimisation to find optimal designs for a 25 kJ laser driver. We perform hydrodynamic scaling on these designs to confirm the achievement of high gain for a 2 MJ laser driver, using 2D simulations including alpha heating effects.

Tables

Table 1. Table 1: Design space used in this work (8 total parameters, 5 laser, 3 target). See equations 1a-f for laser pulse parameterisation.

Variable (units)	Min	Max	Step size
$P_{\mathrm{picket}}$ (TW)	0	25	-
$P_{\mathrm{foot}}$ (TW)	0	12.5	-
$t_{\mathrm{delay}}$ (ps)	50	1000	-
$t_{\mathrm{foot}}$ (ps)	50	1000	-
$t_{\mathrm{rise}}$ (ps)	50	1000	-
$R_{\mathrm{outer}}$ (um)	400	650	1
${\Delta }_{\mathrm{ablator}}$ (um)	2	20	0.5
${\Delta }_{\mathrm{ice}}$ (um)	10	100	0.5

Table 2. Table 2: Performance of the optimal designs at different laser energy scales. Note that full hydrodynamic equivalency re-tuning was not performed when hydro scaled to 2 MJ. Burn off/on refers to simulations excluding/including alpha particle heating.

Scale	1D/2D	Design	${\chi }_{S,no\alpha }$	$Y$	DT n Yield	$⟨\rho$R$⟩$ (g/cm2)	$⟨$T${}_i{}^{}⟩$ (keV)	$v_{imp}$ (km/s)	Bang Time (ns)
25 kJ	1D	$O_{1D}$	2.0	3.52	5.49 $\times$ 1014	0.370	4.35	430	3.67
		$O_{2D}$	1.85	3.40	4.10 $\times$ 1014	0.378	4.16	416	3.50
	2D	$O_{1D}$	1.65	2.02	2.29 $\times$ 1014	0.266	4.22	434	3.67
		$O_{2D}$	1.72	3.05	3.50 $\times$ 1014	0.348	4.08	416	3.50
2 MJ (Burn Off)	2D	$O_{1D}$	-	-	5.86 $\times$ 1016	1.20	4.84	433	15.80
		$O_{2D}$	-	-	1.02 $\times$ 1017	1.46	4.85	416	15.05
2 MJ (Burn On)	2D	$O_{1D}$	-	-	1.02 $\times$ 1018	1.00	7.89	433	15.92
		$O_{2D}$	-	-	2.24 $\times$ 1019	1.16	18.18	416	15.12

Table 3. Table 3: Optimal designs at the 25 kJ scale derived from 1D ( O 1 ​ D O_{1D} ) and 2D ( O 2 ​ D O_{2D} ) simulation based optimisation.

Variable (units)	$O_{1D}$	$O_{2D}$
$P_{\mathrm{picket}}$ (TW)	3.19	4.53
$P_{\mathrm{foot}}$ (TW)	0.65	0.80
$t_{\mathrm{delay}}$ (ps)	806	757
$t_{\mathrm{foot}}$ (ps)	782	668
$t_{\mathrm{rise}}$ (ps)	1000	995
$R_{\mathrm{outer}}$ (um)	436	429
${\Delta }_{\mathrm{ablator}}$ (um)	9.5	9.5
${\Delta }_{\mathrm{ice}}$ (um)	70.5	76.0

Equations

$$P(t)$$

$$P_{\mathrm{picket}}(t)$$

$$P_{\mathrm{main}}(t)$$

$$\tau$$

$$x$$

$$\xi$$

$$\Phi =\frac{\rho R}{6.3{g/cm}^2+\rho R}$$

$$Y$$

$$+\frac{1}{\widehat{\rho R}}\max (\rho R-\rho R_{\chi =1},0),$$

$${\chi }_{S,no\alpha }$$

$$\rho R_{\chi =1}$$

$$\widehat{\rho R}$$

$$S$$

$$\rho R$$

$$M_{\mathrm{stag}}$$

$$⟨P⟩$$

$$P(y_1|y_2=y^{\ast })$$

$${\mu }_1=\rho y^{\ast },$$

$${\sigma }_1=\sqrt{1-{\rho }^2}\sigma ,$$

$$\rho \propto k(x_1,x_2),$$

$$k_{\nu }(\mathbf{x},{\mathbf{x}}^{\prime })=$$

$$\rho =$$

$$\mathrm{LEI}(X)$$

$$Z$$

$$R$$

$$t$$

$$E$$

$$f_{\alpha }$$

$$Y_i<\mu (X_i)+z_{\alpha }\sigma (X_i)],$$

$$z_{\alpha }$$

$$Y_i(X)$$

Full text

I Introduction

Inertial confinement fusion (ICF) seeks to compress and heat deuterium–tritium (DT) fuel to thermonuclear conditions by coupling driver energy into a spherically imploding shell [1, 2]. In laser direct drive[3] (LDD), multiple laser beams directly illuminate a capsule, generating an ablative rocket that accelerates the shell inward. Achieving ignition[4, 5, 6] and appreciable burn propagation requires a careful orchestration of target dimensions and laser pulse history to control shock timing, adiabat, and implosion velocity, while mitigating hydrodynamic instabilities [7, 8]. These competing considerations create a high-dimensional, strongly coupled design space in which small changes to one parameter often necessitate compensating changes in other design parameters.

Radiation–hydrodynamics (rad-hydro) simulations are indispensable for navigating this space, but their computational cost grows steeply with dimensionality. One-dimensional (1D) simulations enable large surveys and provide an upper bound on performance through perfect spherical symmetry, whereas multi-dimensional (2D/3D) simulations are needed to capture Rayleigh–Taylor instability growth that degrades fusion yield[9, 10, 11]. However, the computational cost of 2D calculations can exceed that of 1D by two orders of magnitude, and 3D simulations are more expensive still.

These characteristics motivate multi-fidelity (1D, 2D or 3D) and active-learning machine learning strategies that extract as much information as possible from inexpensive low-dimensional simulations while judiciously deploying costly high-dimensional simulations. Surrogate models can learn mappings from design variables to outcomes and provide predictive uncertainties for Bayesian optimisation [12]. Multi-fidelity extensions exploit correlations across fidelities to accelerate convergence [13, 14]. However, standard multi-fidelity formulations often assume modest computational cost ratios between the fidelities and noise models that are uniform across the fidelities. Thus for multi-fidelity ICF design, there is a need for methods that (i) respect extreme cost disparities, (ii) remain robust to non-uniform, fidelity-dependent noise, and (iii) integrate seamlessly with pre-existing rad-hydro codes; which require High Performance Computing (HPC) clusters.

In this work we develop an automated, simulation-based design framework that couples multi-fidelity surrogate modelling with active learning and Bayesian optimisation for LDD implosions. Our contributions are:

A 1D and 2D radiation-hydrodynamics simulation pipeline allowing straightforward coupling to machine learning (ML) models. For 2D simulations we include fixed (beam-mode) and randomised (ablator surface structure) perturbation sources to reflect shot-to-shot variability. 2. 2.

Multi-fidelity surrogate models that transfer information from large 1D ensembles to data-scarce 2D. These include a neural network ensemble with transfer learning and calibrated predictive uncertainties, alongside Gaussian process baselines. 3. 3.

An active-learning scheme based on probabilistic threshold sampling that concentrates expensive simulations in promising regions of the design space. 4. 4.

Bayesian optimisation at fixed fidelity using multi-fidelity-trained surrogates to find optimal designs in both 1D and 2D. We define optimality via a physics-informed objective that blends an ignition criterion[15, 16], $\chi_{\mathrm{no}\ \alpha}$, with post-ignition burn propagation via areal density[17] $\rho R$, enabling sub-scale optimisation that targets ignition and high gain at a larger hydrodynamic scale.

We demonstrate the framework on an OMEGA-relevant LDD problem with the OMEGA[18] 60-beam geometry, 25kJ driver and an eight-parameter design space (laser picket/foot/rise and layered CH/DT target). We first amass a large 1D database to train surrogates and narrow down the search domain. We then conduct a targeted 2D campaign informed by active learning. We find that transferring structure learned in 1D dramatically reduces the number of 2D simulations required to locate robust designs. Importantly, we show that the 2D-optimised design exhibits substantially improved hydrodynamic stability compared to the 1D-optimised counterpart at the same energy scale. Finally, we assess performance under hydrodynamic scaling[19] to 2MJ (National Ignition Facility[20]-scale) with both burn-off (no alpha heating) and burn-on (with alpha heating) simulations to gauge ignition margin and post-ignition burn.

The remainder of the paper is organised as follows. Section II presents the methodology, defining the design space, objective, and simulation models and detailing the surrogate, active-learning, and optimisation components. Section III presents the optimisation studies in 1D and 2D. Section IV reports the performance of the optimal designs at 25kJ and at 2MJ with burn-off/on analyses. Future work and conclusions are then discussed.

II Methodology

II.1 Design

Design of laser direct drive ICF implosions involves the prescription of a laser pulse (incident laser power as a function of time) and a target (spherical layered target dimensions)[3]. Each of the design parameters (both individually and in combination) are strongly linked to behaviours of the implosion. By tuning of the parameters, one can alter the implosion design (implosion velocity, shell mass, etc.) to achieve the desired conditions at stagnation. However, this is no simple task due to the complex non-linear relationships created by the many physicial processes inherent to ICF. Human designers must consider the qualitative effect of each design decision and use simulations as a tool to quantify these effects. In this work we look to develop a framework that can learn these relationships and thus autonomously move towards an optimal design.

To proceed we must parameterise the laser pulse and target design. The full LDD design space is very high dimensional and therefore we look to reduce the problem down to $<10$ parameters. Firstly, we consider LDD implosions with a 60-beam 25 kJ laser driver. This is relevant to the OMEGA laser facility[18], allowing comparisons between the designs presented in this work with the large OMEGA experimental and simulation database. For the pulse design, we consider a single picket, a foot and power law rise to peak power ($P_{\mathrm{peak}}$ = 25 TW). For the power law rise we use a simplified form of the functional prescription in Gopalaswamy et al.[21]:

[TABLE]

Where $\sigma_{t}=25$ ps and $\beta$ = 9/5. For the target, we use an undoped CH ablator layer, a DT ice layer and a central DT gas region; parameterised as two thicknesses and an outer radius. This laser and target design has a total of 8 free parameters. We must also pick acceptable ranges for these parameters, and for the target parameters we must pick discrete increments such that the initial target fits precisely on an Eulerian grid. Table 1 details the design space, note that all pulse duration parameters were capped at 1 ns to prevent excessively long pulse durations.

With a design space chosen, one must define an objective function. An objective function is used to quantify progress towards some optimal state, where the objective function is maximised. Ultimately, the goal of ICF is to optimise fusion yield, however for sub-scale implosions, which cannot ignite, a different choice of objective is often needed. In the literature, it is common for sub-ignition scale implosions to optimise an ignition metric which can be hydrodynamically scaled[19] to higher laser energies, such as $\chi_{no\ \alpha}$ [15, 16]. Once an ICF target has ignited, then burn propagation into the fuel mass begins. For burn propagation, the metric of success is the areal density ($\rho R$). This can be seen through Fraley’s formula[17] for the fraction of DT fuel, $\Phi$, which undergoes fusion:

[TABLE]

We can therefore aim to design implosions at a lower energy scale which would ignite and achieve high gain at a higher energy scale. The ICF optimisation problem can be broken up into sub-problems to tune aspects of the implosion to achieve the overarching goal [21]. However, we will focus on the global optimisation problem (as in Ref. [22]) and therefore need only provide a single objective function. Guided by the need for ignition and burn propagation at the higher energy scale, we devised an objective which aims to smoothly combine criteria for ignition and burn propagation into a single scalar function:

[TABLE]

Where $S$ is the hydrodynamic scale factor, $Y_{DT}$ is the DT yield, and $M_{\mathrm{stag}}$ is the stagnated mass. The objective function aims to achieve ignition at scale $S$ (occurring when[16] $\chi_{S,no\ \alpha}\sim 1$) before switching focus to increasing areal density above that required for ignition. In this work, we will consider scaling laser driver energies from 25 kJ (OMEGA scale) to 2 MJ (NIF scale). This leads to a scale factor $S\sim 4.3$.

To calculate the quantities required for $\chi_{no\ \alpha}$ and $\rho R_{\chi=1}$, we define the following concrete multi-dimensional definitions for implementation within the radiation-hydrodynamics simulations (in any dimensionality):

[TABLE]

Where $V_{\langle P\rangle e^{-2}}$ is the volume of fuel with pressures exceeding 2 e-foldings of the burn-averaged pressure, $\langle P\rangle$.

II.2 Radiation-hydrodynamics simulation

The design code used in this work was Chimera, a fixed-grid Eulerian radiation magnetohydrodynamics code, details of ICF implosion modelling with Chimera can be found in the literature[23, 24, 11, 25, 26, 27]. To model the laser drive, Chimera was coupled to SOLAS, a 3D laser raytracing model, the details of which will be given in a future publication. Finally, a number of additional code developments were employed to speed up simulation time, namely implicit methods for both flux-limited thermal conduction and multi-group flux-limited radiation diffusion, physics for which Chimera has traditionally used explicit methods[28, 27]. The atomic physics model SpK[27] was used to generate multi-group opacity tables. To facilitate automation, control has been shifted to external input decks and standard decks developed for laser direct drive applications.

The purpose of 2D simulations is to capture hydrodynamic instability growth which degrades performance over the perfect spherical symmetry of 1D. The laser beam geometry introduces a seed for this instability growth in the form of low mode intensity variations on the capsule surface, known as beam mode[29, 30]. Laser speckle operates at much shorter lengthscales than beam mode, introducing high mode sources of instability[31, 32, 33]. Capsule defects, like surface roughness and voids, are also key seeds for instability growth for both low and high mode [34]. To capture some of these effects in our 2D simulations we (a) used the OMEGA beam port geometry to introduce beam mode and (b) introduced density perturbations on the outer 2.5um of the CH ablator. These density perturbations could have been kept constant across different simulations; however, we opted to reseed the random perturbations to introduce noise to the simulated outputs. This gives better correspondence to the reality of shot-to-shot variation in experiment and ensures our methods are robust to such noise. The level of perturbation introduced led to $\sim$1% RMS areal density variations for the initial conditions.

For multi-dimensional Eulerian simulations of spherical implosions it is beneficial to run in spherical geometry for as long as possible to prevent non-spherical grid-based artifacts. However, multi-dimensional regular spherical grids have small cell sizes towards r = 0, requiring small numerical time steps for stability. A restart to other geometries can be used to eliminate this restrictive time step requirement. For this work, an automatic restarting scheme was devised whereby initially 2D spherical ($r-\theta$) geometry was run up until the first shock approached r = 0. At this time, the variables for restart are dumped to disk and a restart is initialised onto a 2D cylindrical ($r-z$) grid. The outer extent of the cylindrical grid was found by identifying the outermost radius at a few percent of critical density.

The above code optimisations allowed run times of approximately 1 and 200 core hours for 1D and 2Ds respectively. For OMEGA scale implosions, a radial resolution of 500nm was used and 256 polar cells in 2D. The simulations were run on 32 and 128 CPU cores for 1D and 2D respectively.

II.3 Surrogate modelling

Surrogate models are a key component of machine learning and optimisation methods. These models use sample data to learn an approximate functional relationship between the inputs ($X$) and the outputs ($Y$). Once trained on this data, the surrogate model is much cheaper to sample from than creating new sample data, i.e. cheaper than running a rad-hydro simulation. Additionally, surrogate models which provide an uncertainty estimate in their prediction are particular useful, allowing Bayesian optimisation [12].

II.3.1 Gaussian Processes

Gaussian processes (GPs) are a common surrogate model for optimisation purposes. There have been a number of works using GPs as surrogate models in ICF modelling and optimisation [21, 35].

For introductory purposes, let us consider a simplified example involving two pairs of variables $x_{1},y_{1}$ and $x_{2},y_{2}$. We expect some relationship between $y$ and $x$, with noise, such that $y=f(x)+\epsilon$. However, we have no functional form for $f(x)$. To proceed, we will assert that the joint distribution of $y_{1}$ and $y_{2}$ is Gaussian with mean 0 and equal variances $\sigma^{2}$. Now let us consider the conditional probability of $P(y_{1}|y_{2})$, given that $y_{2}$ now has some known value, say $y^{*}$:

[TABLE]

where $\rho$ is the correlation coefficient between $y_{1}$ and $y_{2}$. Within this toy model, our conditional mean and variance prediction on $y_{1}$ depends only on the degree of correlation between $y_{1}$ and $y_{2}$. Gaussian Processes make a statement on this degree of correlation by proposing that:

[TABLE]

Where $k(x_{1},x_{2})$ is the kernel function. In other words, GPs say (depending on the choice of kernel) that data points’ proximity in $x$ dictates their level of correlation in $y$. In practice, GPs handle arbitrarily many inputs (often $x\in\mathcal{R}^{d}$) and produce real-valued function outputs. The hyperparameters of the GP kernel, as well as $\sigma$, are found by likelihood maximisation over the observed data, which can be noisy. We direct the interested reader to references which give a full background on GPs [36, 37].

For our purposes, we must make choices on the GP kernel[38] and then fit the GP to the training data. For single-fidelity settings (all the training data comes exclusively 1D or 2D simulations), we used the Matern kernel[37] with anisotropic lengthscales (i.e. a lengthscale per design parameter):

[TABLE]

Where $l_{i}$ is the $i$-th lengthscale, $\nu$ is the smoothness parameter (we use 5/2), and $K_{\nu}(x)$ is the modified Bessel function of the second kind of order $\nu$. Fitting to the training data yields best fit values for the lengthscales $l$ and signal variance $\sigma^{2}$.

However, GPs have some key limitations, a key issue being their poor scaling with training data size. For $\mathcal{O}$(10k) samples, training the GP can become excessively slow. In practice, we found training a GP surrogate with close to the full dataset became slower than performing a Chimera simulation, even when using a GPU for GP training ($\sim$ 10 min training time). This was particularly problematic for multi-fidelity (i.e. using 1D and 2D data) GPs where we modified the kernel to introduce discrete fidelity dependence. In these scenarios, we needed an alternative surrogate model which scaled well with training data volume.

II.3.2 Ensemble of Neural Networks

Neural networks (NNs) are universal function approximators and can therefore act as surrogate models, with previous examples already in the field of ICF[39, 40, 41]. Similar to GPs, NNs use training data to learn a functional relationship between the inputs ($X$) and the outputs ($Y$). In this work, we used a multi-layer perceptron (MLP) NN architecture which involves a sequence of layers of fully-connected “neurons”, the basic unit of a MLP. Neurons apply a linear transformation to the incoming data before applying a non-linear activation function. Designing a good MLP architecture involves the tuning of a modest number of hyperparameters, for example: number of layers, number of neurons per layer, activation function, learning rate, and regularisation parameters.

Unlike Gaussian Processes, NNs do not inherently provide any uncertainty estimate in their prediction of $Y$ given $X$. A common method to model uncertainty with NNs is the ensemble approach [42]. For this, a number of NNs are trained using different random splits of the training data. This causes their predictions to vary, giving a range of results. This is used to quantify the uncertainty in the prediction. To ensure meaningful uncertainty predictions, we performed a calibration scaling on the ensemble predictions, as detailed in Appendix A.

An ensemble size of 25 was used for the 1D and 2D ensemble models. The MLP model was written using PyTorch[43], which allowed its easy integration in the Bayesian optimisation framework BOTorch[44], see Section II.5. A combination of manual and automated tuning (number of layers, neurons, batch size, learning rate) was used to find the final architecture, shown in Fig. 1. Training of each member of the ensemble used the AdamW gradient-based optimiser[45] to minimise the mean squared error between prediction and truth values in the training data. Weight-decay (value of 0.05) was used to regularise the model and prevent over-fitting to the training data. A learning rate schedule (starting at 0.01 and halving every 20 epochs) was found that allowed rapid training in 80-100 epochs (passes over the full training data set) for the 1D surrogate.

Another key advantage of NN surrogate models is the ease of transfer learning[47]. In the context of this work, this means an ensemble of NN surrogates trained on 1D simulation data can be transferred to predict the outcome of 2D simulations, even with limited 2D data. The transfer learning process proceeds as follows:

Train an ensemble of NN surrogates on randomised train-test-validation (60:20:20) data splits of the 1D data set 2. 2.

The first N layers of the NN surrogates are ‘frozen’, such that their trained parameters will not change in subsequent training 3. 3.

Perform transfer learning by re-training the ensemble of NN surrogates on randomised train-test (70:30) data splits of the 2D data set

In this manner, information from the 1D (correlations between input parameters, useful non-linear embeddings of the inputs, etc.) will be transferred to a 2D surrogate. In practice, it was found the freezing all but the last 2 layers produced good results when tested on validation data.

II.4 Active learning

The input parameters for laser direct drive (target dimensions, laser pulse parameters) are highly correlated in good designs. For example, shocks launched by the laser are optimal if they can be merged together precisely[48], this correlates the thicknesses of target layers to laser parameters like picket powers and durations. This correlation leads to good designs residing within small, complex volumes in the design/input space. To enable effective optimisation, we wish to populate these volumes with more sample data rather than uniformly filling the input space[21]. We opted for an active learning[49] approach to achieve this.

For active learning, we will use a surrogate model trained on uniformly sampled data to narrow down the search space for good designs. Based on the surrogate model predictions, new sample data will be taken (via rad-hydro sims) and used to update the surrogate model. In this way, we can increase the sample frequency within regions we care about, while training the surrogate simultaneously. This also has the additional benefit that batches of candidate points with high objective function values can be pulled from the surrogate model and run concurrently, i.e. in batches. In this work, we employed the active learning algorithm described in Algorithm 1. This shares features with canonical level set estimation[50] algorithms.

Where $\Phi\left(X\right)$ is the cumulative normal distribution and $\mathcal{U}(a,b)$ is the random uniform distribution between $a$ and $b$. There are a number of hyper-parameters to be chosen, $P_{\mathrm{min}}$ and $\hat{Y}_{\mathrm{threshold}}$. It was found that $P_{\mathrm{min}}$ values between 0.25 and 0.5 gave a good balance between exploration and exploitation. The threshold is chosen such that designs with expected outcomes better than the threshold are searched for.

II.5 Bayesian optimisation

Initial sampling and active learning was used to build a comprehensive data set to train a sufficiently accurate surrogate model to be used in Bayesian optimisation. Bayesian optimisation defines an optimisation problem in terms of maximising an acquisition function, which itself is a function of the surrogate model mean and uncertainty predictions, $\mu(X)$ and $\sigma(X)$ respectively. This optimisation problem over the acquisition function is computationally cheaper than the ‘true’ optimisation problem over the simulator outputs.

We used the same Bayesian optimisation strategy throughout this work. Firstly, we restrict to simulation batch sizes of 1 i.e. sequential optimisation. This allows for the most well defined acquisition function and optimisation problem. Secondly, we used the log Expected Improvement[51] acquisition function, which is defined as:

[TABLE]

Where $\Phi\left(X\right)$ and $\phi\left(X\right)$ are the cumulative and probability distribution functions of the standard normal distribution. This acquisition function is maximised w.r.t. $X$ using gradient-based methods. A new simulation is then ran at the optimal $X$ found. This produces a new $X,Y$ pair for the dataset, the surrogate is retrained and the process is repeated. The numerical solution to this optimisation problem was performed by the BOTorch library[44] and the coupling to the simulator is handled by software described in the following section.

It should be noted that we have not taken a fully multi-fidelity Bayesian optimisation approach[52], even if we are using multi-fidelity surrogates. This is because of the large difference in computational cost ($\sim$ 200x) between our fidelities (1D or 2D rad-hydro sims). In traditional multi-fidelity Bayesian optimisation, the optimiser chooses the next candidate $X$ as well as the fidelity based on the information gain per unit cost. Due to our cost disparity, 2Ds are run very sparingly so cost weighting would need to be used. Instead, we opted to run Bayesian optimisation at a fixed fidelity but using surrogates trained on both fidelities. This allows the transfer of information between fidelities via the surrogate but updates to the database are strictly within a single fidelity.

II.6 Orchestration

An orchestration framework is required to efficiently couple together the machine learning and simulator models described above. The job of orchestration software is to: convert candidate sample points into valid simulator inputs; execute and wait for batches of simulations to run within a HPC environment; post-process the simulations for valid training data; re-train the surrogate model and finally query the explorer/optimiser for new sample points. Fig. 2 summarises this workflow.

While the specifics of this workflow are simulator and surrogate model dependent, the software package developed for this purpose, mille-feuille, is open-source[53] and aims to abstract this workflow process.

Creating automated design programs in this way allows multiple levels of parallelism to be exploited. At the lowest level, the simulator is parallelised by MPI with load balanced domain decomposition. Within the automated design process, batches of the simulator can be launched simultaneously to use all the allocated HPC resource. At the top level, multiple instances of the automated design program can be launched as an array, learning from the same centralised data sets.

III Optimisation Study

III.1 1D design

An extensive initial sampling 1D simulation ensemble was run due to the relatively low cost of 1D simulations. The input domain defined in Table 1 was transformed to the unit-hypercube and the quasi-random Sobol sequence[54] was used as a space filling design[36]. Of order 12.5k 1D simulations were ran based on this Sobol design. As expected the vast majority (97%) of these simulations represented poor designs (quantified as objective $Y<1$) which are unlikely to ignite at the 2MJ scale. However, of those simulations with $Y>1$ the range of the input parameters still covered the whole design space. This is because the inputs for this dataset of “good” designs are correlated. For example, ice and ablator thickness are anti-correlated with each other (-0.6 Pearson coefficient) and outer radius (-0.16 and -0.35 respectively), suggesting good designs lie within a certain target mass range. For the laser pulse, picket power is positively correlated with foot power (+0.35) and negatively with the delay between picket and foot (-0.3). This reflects the importance of the front end of the pulse in shock timing [48]. The distribution of objective values from this initial sampling can be seen in Fig. 3.

The comprehensive initial sampling allowed the training of a single fidelity GP surrogate for active learning. Since only 3% of the initial sampling lay within the region of interest ($Y>1$), active learning with a threshold of 1 and $P_{min}$ = 0.25 was used to find $\sim$ 1k points within this region of interest. As expected, the active learning sample shows similar statistical correlations as the thresholded initial Sobol samples. The active learning step increased the population of high objective function observations, increasing surrogate model confidence in this region.

Finally, Bayesian optimisation was performed to find the optimal design in 1D within the provided design space. Two approaches were taken. Firstly, a GP surrogate model was used, however training time was becoming obstructive due to the large 1D dataset. For the GP, the training data was downselected by reducing the design space to encompass all points with $Y>1.5$, this reduced the number of training samples from $\sim$ 14k to $\sim$ 5k. Secondly, the NN ensemble was used, for which the whole 1D dataset was kept. Both surrogates were used to optimise Log Expected Improvement. All but the last two layers of the NNs were frozen for retraining during the optimisation loop. The optimisations with the different surrogates converged towards similar optimal designs ($Y$ values within 0.1 of each other), however it was the NN surrogate which found the highest performing design after $\mathcal{O}(10)$ new samples.

III.2 2D design

It is evident that any target optimisation using only 1D simulations will produce optimistic designs, as the designs do not need to be robust to hydrodynamic instabilities to achieve high performance in 1D. However, there are great benefits to transferring some of the information learned from 1D simulation ensembles to optimisation using multi-dimensional radiation hydrodynamics simulation. In particular, 1D results represent the upper limit of performance for 2D calculations. Therefore, by transferring information about the design space from 1D to 2D we can reduce the search space drastically.

Therefore, our 2D design process proceeded as follows. A subset of the inputs from 1D active learning samples were re-run as 2D simulations, we will refer to this as our “transfer sample”. This transfer sample allowed the training of 2D surrogate models. While both multi-fidelity GPs and NN ensembles were tested, it was found that the NN ensembles performed equally well with considerably lower training times. Trained on the 1D dataset and the 2D transfer samples, surrogate model predictions for both 1D and 2D were compared for inputs ($X$ values) in the 1D dataset. In Fig. 4, we can see the 2D surrogate has learnt from the transfer sample that performance is degraded relative to 1D. It is also observed that the degradation is a strong function of input parameters, particularly ice thickness for these simulations. It is worth repeating that the 2D data also is inherently noisy as different random surface perturbations are applied for each simulation. This is reflected in the large surrogate uncertainty shown in the error bars, noting that these uncertainties have been calibrated using the method shown in Appendix A.

Approximately 128 2D simulations were run within an active learning loop to explore the space of good 2D designs. For the active learning, a threshold objective value of $2$ and $P_{min}$ of 0.5 were used, aiming to exploit the space rather than to explore. Fig. 5 shows the new samples created by the active learning method. This clearly shows the transfer learning has been effective in narrowing down the search space for 2D designs based on the 1D dataset. The 2D surrogate directs samples towards large objective values and shows no examples of poor design choices (lowest objective value $\sim 0.9$). The transfer sample does not contain poor 1D designs, $Y<1$, so it is clear that the understanding of the design space found in 1D is maintained in the transfer learning, in particular what regions of design space to avoid due to low performance.

Finally, Bayesian optimisation was performed using the 2D NN ensemble surrogate only, with transfer learning performed using all available 1D and 2D data. With a small number of runs, the algorithm found a small region of input space with robust 2D performance, with a maximal objective function of approximately $3.0$, to be compared to $3.5$ in 1D.

IV Results

In the following two sections we will provide a detailed examination of the performance of the optimal designs at both the design scale (25 kJ) and scaled up (2 MJ). Table 2 summarises the key performance parameters for an overview of the design performance.

IV.1 Performance of optimal designs at 25 kJ scale

We will first consider the result of the 1D and 2D design optimisation at the design energy scale of 25 kJ. We will assess the design performance using the objective function, Eq. 3a, as well as key ICF performance parameters[55] (DT fusion yield, ion temperature, implosion velocity, etc.). The optimal design parameters are given in Table 3 and the key performance parameters are summarised in Table 2.

The relative differences between the 1D and 2D optimal designs (denoted $O_{1D}$ and $O_{2D}$ respectively) give rise to only small differences in performance in 1D, with areal densities within 2% of each other, with the $O_{1D}$ design having a higher $\chi_{S,\mathrm{no}\alpha}$ value. To understand these differences in the designs we will look at important implosion characteristics related to both 1D performance and hydrodynamic stability. By inspecting the pressure gradient shown in Fig. 6,

we can see the optimiser achieves a reasonable degree of shock timing for both designs. It is worth noting that this is not specifically optimised for, but arises when maximising our chosen objective. This allows these designs to achieve high areal densities (370 - 380 mg/cm2). These areal densities are considerably higher than has been achieved in experiment at OMEGA[56] ($\lesssim$ 250 mg/cm2), pointing to the optimism of 1D design calculations.

From the 1D simulations, we can also look at key metrics of hydrodynamic stability. These are the inflight aspect ratio[57] (or IFAR) and the minimum adiabat[2] of the DT fuel, $\alpha_{DT}$. These stability metrics are plotted as a function of time in Fig. 7.

Higher IFAR means the imploding shell is thin in-flight and therefore more easily disrupted by the non-linear growth of Rayleigh-Taylor. The $O_{1D}$ and $O_{2D}$ have peak IFARs of 37 and 30 respectively. Another key stability metric is the shell adiabat (defined the ratio of the total pressure to the degeneracy pressure). Low adiabat shells are more compressible leading to higher final densities at stagnation. However, low adiabats are unfavourable for hydrodynamic stability[58, 59]. Fig. 7 shows that $O_{2D}$ has a consistently higher adiabat than $O_{1D}$. This is achieved by using a higher picket power (4.5 TW vs 3.2 TW), as well as a degree of shock mistiming seen in Fig. 6.

As expected from analysing the stability in 1D, multi-dimensional Chimera simulations show that the $O_{1D}$ design degrades far more than the $O_{2D}$ design in the presence of hydrodynamic instabilities. This is particularly evident in the bangtime density and temperature distribution as shown in Fig. 8. The DT shell in the $O_{1D}$ design is destroyed by the growth of hydrodynamic instabilities, whereas the $O_{2D}$ design keeps the shell relatively intact. The OMEGA beam mode (averaged into 2D $r,\theta$ geometry) is particularly evident in the $O_{1D}$ fuel morphology. It is clear the $O_{2D}$ design (while similar in many design parameters to $O_{1D}$) is more robust to hydrodynamic instabilities.

The trained 1D and 2D surrogate models are also a useful tool to visualise and understand the effect of a single design decision. We can optimise over the surrogate model’s (mean) prediction of the objective, while holding one input parameter at a fixed value. The results from this investigation are shown in Fig. 9. These results show some intuitive findings when comparing the predictions for 1D and 2D. For example, in 2D decreasing picket power below a threshold rapidly decreases performance. Low picket power creates a low adiabat implosion which is more prone to hydrodynamic instability growth. In a similar vein, thinner ice is more disrupted by the growth of ablator surface perturbations - the key source of degradation in our 2D simulations. Therefore, we can see from these trends that the 2D surrogate has found regions of input space which are robust to the level hydrodynamic instabilities that was seeded. Amplifying the level of ablator density perturbation applied would push the trends further, heading towards smaller capsules, with thicker ice and higher picket power/adiabat.

IV.2 Performance of optimal designs at 2 MJ scale

In this work, the design objective was devised such that the optimal targets hydrodynamically scaled to a 2 MJ laser driver would ignite and produce high gain. In particular, via a multi-fidelity approach, we aimed to find targets that were robust to hydrodynamic instabilities seeded at the ablator surface. In this section, we investigate the scaled design performance via direct simulation at the 2 MJ scale with alpha heating physics included. We used the particle Monte Carlo burn package in Chimera as described in Tong et al.[25]. This model was extended to include the additional effect of fuel depletion via a burn fraction advected around with the DT material which was used to reduce the reaction rate.

Hydrodynamic scaling by scale factor $S$ involves the following dimensional scalings[19]:

[TABLE]

Where $R$, $t$ and $E$ are the spatial, temporal and energetic scales. However, due to non-hydrodynamic physics (e.g. radiation and thermal transport) naively applying this scaling to initial conditions does not lead to hydrodynamically equivalent implosions. Nora et al.[19] note additional changes to the laser pulse and target dimensions are required to restore hydrodynamic equivalency, mainly exchange of ablator mass for DT ice. Here we will not perform this re-optimisation but instead apply hydro-scaling to the 2D simulations at a restart time approaching stagnation, making the assumption that a hydro-equivalent implosion can be achieved up to this time.

First, both designs have sufficient ignition margin to ignite in 2D at the 2 MJ energy scale. However, the yield amplification due to alpha heating (defined as the ratio of neutron yields burn-on to burn-off) are very different with values of 17 and 217 for the $O_{1D}$ and $O_{2D}$ designs respectively. From the burn-off areal densities, one might expect ignited designs to achieve burn fractions of 16% and 19% respectively, based on Eq. 2. However in the burn-on simulations, burn fractions of 0.9% and 15% were achieved. This shows canonical 1D burn propagation was not achieved in the $O_{1D}$ design due to large asymmetries in the confining shell.

Fig. 10 shows the temporal evolution of the burn characteristics for the two designs. Comparing burn-off and -on results for the same design shows higher burn rates and temperatures but lower areal densities. The lower areal densities are due to the rapid expansion driven by increased hotspot pressure from alpha heating. Now contrasting the two designs, the strongly burning $O_{2D}$ design shows super-exponential growth the burn rate, shown by the solid red curve which has positive second derivative in the rise on the log-linear scale. Burn-averaged ion temperature also grows rapidly during this time up to a peak of 22 keV. During the fall of the burn ($t>15.1$ ns), burn-averaged ion temperature shows a second increase as shock heating of the remaining low density, free-falling shell[60] produces large ion temperatures (10s keV). This resembles the hohlraum re-heating[61] phenomena observed on the National Ignition Facility for ignition experiments. However this free-falling DT heating does not contribute significantly to increased burn. The $O_{1D}$ burn-on simulation shows less pronounced increases in ion temperature and burn rate reflecting the truncated burn propagation. This can be attributed to poor hotspot confinement from thin regions in the shell. This allows the hotspot to form aneurysms through holes in the shell[62, 63].

The differing degree of burn propagation leads to very different conditions at the time of peak neutron production, as shown in Fig. 11. The strongly burning $O_{2D}$ design has become a near uniform ball of burning and expanding DT. High localised burn up fraction, $\Phi\sim 20\%$, leads to reduced reaction rate, $R_{DT}$, in the core of the burning plasma, where temperatures exceed 30 keV. In contrast, the $O_{1D}$ design shows a largely perturbed shell into which burn has not propagated. Hotspot fuel has blown out through holes in the shell showing it has lost confinement before burn of the dense fuel could be established. This results in the $O_{1D}$ design having far lower performance than the $O_{2D}$ design, as discussed above.

It is clear from these direct numerical simulations at the 2 MJ scale that the objective function used (Eq. 3a) and the multi-fidelity methodology was able to design high gain designs at a higher energy scale than simulated. The optimised 2D design allowed for a more robust burn due to its better stability leading to a more uniformly confining shell.

V Future work

In this paper we describe a framework for performing multi-fidelity optimisation around the Chimera simulation code. This provides many avenues for future work, using Chimera’s full capabilities as a radiation magnetohydrodynamics code. For example, additional fidelities can be introduced such as 3D simulation[23], the modelling of cross-beam energy transfer[64] or the use of experimental data. It is also important to introduce further sources of perturbations, such as DT ice roughness[65, 66] or capsule stalk[67], to ensure optimised designs are robust to instabilities derived from experimentally realistic perturbations. There are also magnetically driven inertial fusion schemes (such as MagLIF) for which optimal designs at higher driver energies are an active area of research[68, 69].

We have operated with a relatively small design space of 8 parameters. Scaling to larger design spaces is necessary to explore a greater variety of possible designs. Some key additional parameters being considered for next generation inertial fusion designs relate to target layers (dopants, foams)[70, 71] and laser parameters (beam size, zooming, bandwidth)[72]. However, the curse of dimensionality means that larger design spaces require larger numbers of sample points to train the surrogate. Additionally, the global optimisation of the acquisition function becomes increasingly difficult in higher dimensions. This will require new optimisation strategies, such as the splitting into smaller sub-problems[21] or using gradient-based optimisation methods directly from the simulator[73, 74].

VI Conclusions

In this work we show that machine learning models can use simulation data from multiple fidelities to optimise designs in the context of laser direct drive inertial fusion. In particular, we show that a large number of 1D simulations, with inherently no hydrodynamic instability growth, could be used to narrow down the search space for 2D simulations. It was found that an ensemble of neural network based surrogate models was effective at transferring this information between fidelities. This surrogate model also enabled active learning and Bayesian optimisation methods to automatically tune the laser and target design parameters. These methods were robust to noise introduced to the 2D simulation results via randomised ablator density perturbations.

The optimal 2D design was more hydrodynamically stable than the optimal 1D design and showed only weak degradation in performance at the 25 kJ energy scale. Hydrodynamically scaled simulations for a 2 MJ laser driver including alpha heating physics showed large differences in the achieved yield amplification, 17 and 217 for the $O_{1D}$ and $O_{2D}$ designs respectively. The $O_{1D}$ design failed to propagate burn into the dense fuel due to loss of confinement from hydrodynamic instability growth. In contrast, $O_{2D}$ design was highly performant at the 2 MJ scale. This demonstrates that our devised objective function and multi-fidelity design framework correctly identified a high-gain scaled-up design, robust to the presence of hydrodynamic instabilities.

Acknowledgements

This research received support through Schmidt Sciences, LLC, the International Atomic Energy Authority (IAEA) AI for Fusion Coordinated Research Project and the Imperial College Research Fellowship program.

This work made extensive use of the Imperial College London RCS HPC systems and the ARCHER2[75] UK National Supercomputing Service (https://www.archer2.ac.uk).

The authors have benefited from many useful discussions with individuals at I-X, the Australian National University and Data61. We also thank Dr Brian D Appelbe for his helpful review of this work.

Appendix A Ensemble Model Uncertainty Calibration

While an ensemble of ML models can give a distribution of predictions at a point, it is not guaranteed that the statistics of these predictions is representative of the data it was trained/tested on. Often the minimisation of a mean squared error loss ensures the mean prediction is close to the data, however the ensemble variance might be miscalibrated.

Following Ledda et al.[76], we aimed to test our ensemble model calibration by computing the frequency at which test data lies within a given confidence interval of the ensemble model. This is defined mathematically as:

[TABLE]

Where $f_{\alpha}$ is the frequency of test data $D$ in confidence interval $\alpha$, $X_{i},Y_{i}$ are data from test data $D$ which is size $N_{D}$, $\mu(X)$ and $\sigma(X)$ are the ensemble model predictions of mean and standard deviation at input $X$, $I$ is the indicator function, and $\Phi^{-1}$ is the inverse cumulative probability function of the normal distribution.

Performing this analysis with both our 1D and 2D surrogate ensembles revealed that the trained models were typically ‘over-confident’ in their predictions, i.e. the error bars of the ensemble was too small. Fig. 12 shows this as the uncalibrated curves lie below the perfect calibration line.

It was found that applying a constant scale factor ($\sim$ 3.5) on the model uncertainties regained close to perfect calibration. This calibration constant ($c_{v}$) was then included into our ensemble model predictions as follows:

[TABLE]

Where $Y_{i}(X)$ is the $i$-th models prediction at input $X$ and $\bar{Y}(X)$ is the mean prediction over the ensemble. This transformation ensures that the model mean is unaffected but the variances are scaled by $c_{v}$.