1 Introduction

Today, there remain three commercial nuclear power plants (four reactors) in operation in Switzerland. Along with the fifth reactor that was shutdown in 2019, a certain amount of spent fuel assemblies discharged from the reactors, hereinafter referred to as Used Nuclear Fuel (UNF), has accumulated over the years. A typical fuel assembly for a Pressurized Water Reactor (PWR) is about 4 to 5 m high, 20\(\times\)20 cm large and weights about 400 to 500 kg; its content in radioactive isotopes varies depending on irradiation time (e.g., fission products such as \(^{99}\)Tc, or actinides such as \(^{239}\)Pu or \(^{241}\)Am). Consequently, after being used in a reactor core for several years, the UNFs need to be cooled, stored and eventually disposed of. Facilities used to this aim, referred to as Spent Fuel Management Systems (SFMS), include wet storage pools, dry transport/storage casks and disposal canisters. In Switzerland, UNF recycling in the form of reprocessed fuel such as mixed oxide (MOX) fuel is no longer applied; instead, the concept of so-called permanent direct disposal into a deep geological repository was adopted through a national sectoral plan approved in 2008 by the Swiss Federal Council. The current plan is for the long-term (virtually infinite) underground storage of more than 12,000 UNFs [1]. For illustration, a model of a canister loaded with four UNFs and placed in a underground tunnel (a few hundred meters deep) is presented Fig. 1. One can see four UNFs inside the steel container, surrounded by the bentonite blocks and backfilling.

Fig. 1
figure 1

Tunnel model for the disposal of used nuclear fuels in a specific canister \(^{\copyright }\hbox {LMS}/\hbox {EPFL}\), Underground laboratory in Mont Terri [2, 3]

It can easily be understood that the disposal of such canisters will be subject to very strict safety and environmental requirements. A key aspect will also be to ensure that all the UNF assemblies can be disposed of without resulting in an unnecessary large number of canisters. To achieve this, an optimization of the canister loadings will therefore be needed, and the scope of this paper is to present a methodology that could broaden the options for such optimization. More specifically, for a given initial population of UNF assemblies, the present optimization will aim at finding loading combinations that allow using the exact number of canisters required to dispose of all the available UNF population, while ensuring strict compliance with the two selected safety and design related constraints, namely the neutron multiplication (“k\(_{\mathrm{eff}}\)”) and the thermal heat output (or “decay heat”) [4,5,6].

To this aim, the principles of the proposed methodology are to take benefit of more knowledge on the UNF isotopic contents through higher resolution neutronic simulations and to combine this with artificial intelligence techniques in order to search the state-space of allowed loadings. The existing investigated methodology is based on the canister loading with homogeneous UNFs (i.e., four UNFs with similar neutronic characteristics, such as the initial fissile material enrichment, the burnup (BU) achieved at the end of irradiation, and their cooling time). The drawback of such approach is that it does not take advantage of possible combinations of different assemblies (for instance with various burnup values), consequently leading to an increase in the number of disposal canisters. Alternatively, the methodology presented in this paper allows to take into account heterogeneous UNFs in the same canister. The effective neutron multiplication factor k\(_{\mathrm{eff}}\) for an heterogeneous canister loading can theoretically be computed using a Monte Carlo transport and depletion code such as SERPENT; however, it requires up to several calculation hours for a single canister loading calculation. As a practical solution, artificial neural networks (ANNs) only require a fraction of a second to evaluate the k\(_{\mathrm{eff}}\) of a specific canister loading, once it has been trained, making it possible to investigate optimization algorithm using heterogeneous loading.

This paper is therefore structured as follows. To start, an overview is provided on the context of the presented optimization along with the associated assumptions and strategy. The next section provides details on the geometrical problem, the various input parameters, simulation tools and the considered approximations; the third section gives details on the ANN and its performances; the fourth section presents the use of the ANN with the genetic algorithm together with the main results of this study, and finally in the last section, we discuss the advantages and limitations of the present approach and the possible improvements to be done for later studies.

2 Context of optimization assumptions and strategy

For all stages of the post-operation lifecycle of the discharged fuel, the management and handling of the UNF in SFMS can be summarized as a very large optimization problem (see, e.g., Ref. [7]). The whole optimization combines a multi-dimensional range of interdependent objectives and constraints: to achieve specific design and technological characteristics, to provide sufficient capacity without increasing the fuel cycle costs, to allow for cost effective SFMS maintenance and operation, to ensure strict compliance with all safety criteria and regulatory requirements, and to fulfil all constraints with regards to environmental protection.

In general, this translates into five main SFMS optimization areas: (1) materials and structural design, (2) thermal performance, (3) nuclear criticality prevention, (4) minimization of radiation emissions and dose rates, and (5) operational and maintenance flexibility including logistics. For the latter, not all UNFs will be available at once and this will thus require proper scheduling and coordination, for instance regarding transport, encapsulation and transfer to the final repository.

Each of these components might moreover comprise sub-areas of optimization. For instance, with regards to nuclear sub-criticality, which refers to the required safety margins to prevent the occurrence of a self-sustained nuclear chain reaction, Burnup Credit (BUC) can upon regulatory acceptance be implemented for the optimization of SFMS loadings. This translates into allowing to take into account for the fuel irradiation during reactor operation, when performing criticality safety evaluations. The reason is that BUC will allow to significantly reduce the amount of fissionable materials compared to calculations based on assuming fresh fuel only.

For the work presented here and dealing with disposal canisters, the optimization will focus only on two of the above optimization components, namely criticality and thermal performance, with the overall goal to search for loading combinations that will allow (1) using the lowest number of canisters required for a given UNF population, and (2) having the decay heat distribution for all canisters as uniform as possible. To this aim, it will be assumed that BUC is allowed.

For criticality, the objective function will be the neutron multiplication factor k\(_{\mathrm{eff}}\) of a canister containing UNFs. A general definition of k\(_{\mathrm{eff}}\) is explained in the following. As some fissile materials are present in the UNF after irradiation (mainly \(^{239,241}\)Pu and \(^{235}\)U), there are possibilities to obtain a system where the nuclear chain reaction is self-sustaining, which is to be avoided (see for instance Ref. [8] for details). Such risk depends on the content of the UNF, their location in the canister, the possible water infiltration over time, and on the canister geometry and materials.

As the content of the UNF varies over time (due to the radioactive decay of some materials), the k\(_{\mathrm{eff}}\) needs to be calculated for the full storage time. In neutronics, the calculated quantity called k\(_{\mathrm{eff}}\) is used to represent the criticality state of a system: if k\(_{\mathrm{eff}}=1\), every fission produces on average one more fission, leading to a constant fission level (power plants operate at k\(_{\mathrm{eff}}=1\)); if k\(_{\mathrm{eff}}<1\), the system cannot sustain a chain reaction, and if k\(_{\mathrm{eff}}>1\), the number of fission reactions increases exponentially. In the present study, an upper sub-criticality limit (USL) is considered: k\(_{\mathrm{eff}}\) must be lower than 0.95, allowing for a safety margin of 5 % below 1, or also mentioned as 5000 pcm (per cent mille) below 1 (this value is generally used in criticality safety studies [9]). It must however be noted that the USL used in real licensing applications will usually be refined in order to take into account burnup (see, e.g., Ref. [6]) as well as uncertainties and/or biases due to the employed methods, data as well as SFMS/UNF structures, geometries and materials.

Regarding thermal performance, the total decay heat (DH), representing the canister heat load and thermal output, will be considered here as second objective function. The DH corresponds to the amount of energy still released by the used fuel, once it is unloaded from a reactor core. It originates from the decay of radioactive isotopes and can occur as long as all unstable isotopes have not decayed to stable ones. If the amount of DH is too high, the canister outer surface temperature might reach values high enough to alter the physicochemical properties and characteristics of the multi-barrier system surrounding the canisters (e.g., bentonite buffer and host rock). In the long-term, this might lead to an accelerated degradation of the geological and engineering barriers and/or create migration paths for radionuclide transport to the biosphere in case of canister wall breaches and/or cracks. Currently, a total of 1.5 kW is being considered in Switzerland [10] as design limit for the canister decay heat; this will therefore be used here as second optimization constraint. This DH constraint will consequently inuence the type of UNF loaded together in the same canister, as the UNF decay heat depends on a few characteristics: initial content of ssile materials (only \(^{235}\)U is considered here), how long it was irradiated in a reactor core (or burnup, quantied in MWd/kgU), and how long it stayed on-site after its last irradiation.

For the study presented here, the k\(_{\mathrm{eff}}\) and DH optimization will be performed at the start of the disposal period after an assumed post-discharge cooling and dry storage time of 30 years. Although longer storage times are actually anticipated in many countries, referred to as “an extended intermediate storage” and which could span between 40 to even 100 years, 30 years was considered as sufficiently representative for the purpose of the optimization study. Moreover, it must be underlined that for certain configurations, the k\(_{\mathrm{eff}}\) could increase during the disposal phase (before decreasing again). Such time-dependent effects on the optimization will not be considered here and sub-criticality will thus only be evaluated after 30 years of cooling, i.e. at start of disposal.

If the total decay heat for a loaded canister can be obtained by summing the individual values of each UNF, the k\(_{\mathrm{eff}}\) is not an additive parameter and it also depends on the positions and orientations of each UNF. Additionally, as explained earlier, one k\(_{\mathrm{eff}}\) calculation is usually performed using a Monte Carlo neutron transport code, which is relatively time consuming (several hours on a normal computer). As a consequence, choosing combinations of UNFs and their orientations to minimize k\(_{\mathrm{eff}}\) require a massive amount of simulations. To illustrate the dimension of the possibilities, a simple case of one canister and 4 UNFs leads to more than 10\(^{4}\) combinations (presuming that one UNF can be radially and axially rotated). For 2 canisters and 8 UNFs, the number of combination is over 10\(^{9}\). It is easily conceivable that with a few hours of calculation per case, loading optimization becomes very challenging, if not impossible.

In order to replace Monte Carlo simulations and to effectively speed up the k\(_{\mathrm{eff}}\) calculation, a convenient solution is to use an artificial neural network (ANN). The advantage of such surrogate method is that once the ANN is trained, then a k\(_{\mathrm{eff}}\) can be obtained only in a fraction of a second, therefore opening the possibility of multiple calculations for the optimization purpose.

This is the first goal of the work described in this paper: how to define a ANN for k\(_{\mathrm{eff}}\) calculations and which performances are then obtained (in terms of bias or maximum deviations)? The second goal pursued here is then the optimization of the canister loading with a genetic algorithm, using an ANN for the calculations of k\(_{\mathrm{eff}}\).

Finally, it must be noted that the optimization of SFMS has been subject of a wide range of previous studies. This also includes canister loading optimization for which combinatorial approaches as well as heuristic methods were for instance explored (see Refs. [11, 12]). Usage of optimization techniques such as genetic algorithms and simulated annealing were also investigated, for instance in the context of minimizing costs related to on-site UNF storage after reactor shutdown [13]. Concerning neural networks, they have been successfully used in the past for regression problems in multiple domains. These range from predicting river daily flow [14] to asphalt concrete stability [15] and include deriving electronic temperature in nuclear fusion experiments [16] and many more. ANNs can be in particular a very powerful tool for regression problems with highly nonlinear relation [15, 17].

In relation to nuclear science and engineering, ANNs and artificial intelligence methods have been applied for loading optimization in relation to in-core fuel management (e.g., Refs. [18,19,20,21]).

The work presented here constitutes thereby a continuation of some of the previous studies but with the particularity that here, the intention will be to take benefit of higher resolution methods for spent fuel characterization and combine these with a dual optimization scheme comprising both an ANN and a genetic algorithm (GA) to perform canister loading optimization.

3 UNF, canister and neutronic simulations

In this section, general descriptions of the UNFs, the origins of their contents and their characterizations are provided. Such models and the calculated k\(_{\mathrm{eff}}\) will be the learning database for the ANN, and in total, about 46,000 different canister loadings are considered.

3.1 Description of the assemblies

For the present study, 212 representative UNFs operated in a Swiss PWR are considered, spanning over the range of initial fissile content enrichments (from 1.9  to 4.39%) and burnup values (from 8 MWd/kgU to 66 MWd/kgU). In this context, only uranium fuel (UO\(_2\)) assemblies are here considered, meaning that other operated fuel designs such as MOX are for the time being not included in the optimization.

The selected enrichments and burnup values are representative of realistic cases, either for UNF being fully burned, or for UNF expected from the last reactor cycle before final shutdown (not fully burned). These UNFs originate from a specific Swiss PWR, and can potentially be loaded into a minimum of 53 canisters. If not all DH and k\(_{\mathrm{eff}}\) values are below the limits, more canisters will have to be used, unfortunately increasing the cost of the long-term storage.

This number of 212 UNFs is sensibly smaller than the expected 12,000 UNFs for final repository [5] (when considering all Swiss power plants, and assuming a base scenario of 60 operating years for all the Swiss reactors). In practice, as power plants are still in operation, a number of assumptions would be necessary to cover the full 12,000 UNF cases. But as the range of UNFs selected for this study is large, it still allows to study the applicability of a ANN. The assemblies considered in this study are all geometrically similar: made of \(15\times 15\) vertical rods (225 in total): 20 of them are empty for control rode guide, and 205 are filled with fuel. These latter ones are initially made from uranium oxide UO\(_2\) with a specific concentration of \(^{235}\)U (or enrichment). All rods have a zirconium cladding.

3.2 Simulation codes for reactor physics

A series of simulation codes are used to calculate the decay heat and k\(_{\mathrm{eff}}\) for loaded canisters, as a function of the selected UNFs and their location. Details are provided in the next sections, the goal being to replace such codes for criticality calculations by a neural network.

3.2.1 Reactor core and high resolution UNF simulations

As mentioned, the selected 212 UNFs were irradiated in a Swiss PWR and their irradiation history as well as their physical characteristics are thus known. The time spent in a reactor core, the cycles of irradiation, as well as other information on the reactor are extracted from the Paul Scherrer Institute (PSI) in-house Core Management SYStem (CMSYS) computational platform for reactor core simulations [22, 23]. Within CMSYS, reference assembly/core models based on the Studsvik codes CASMO [24] and SIMULATE [25] are being continuously developed and validated for all the Swiss reactors and operated cycles. Recently, CMSYS was enlarged to couple the core models to the SNF inventory code [26] also belonging to the Studsvik suite.

On this basis, a computational scheme was established such as to allow reconstructing the isotope inventories and associated post-discharge source terms for each axial layer (typically 30 to 50 vertical segments), for each individual fuel rod in each of the operate fuel assemblies, taking into account the detailed irradiation environment that surrounded the assembly during operation. This constitutes a significant advance in terms of spent fuel characterization, and is defined here as high resolution predictions of realistic UNF contents. In addition, this opens perspectives for new optimisation alternatives.

To illustrate this, Fig. 2 presents two simplified views of the central cross section of canisters loaded with 4 different UNFs. The UNFs were selected arbitrary and for demonstration purposes only. The same assemblies are inside both canisters, only their radial rotations are different. Numbers indicate the initial enrichments, as well as their burnup values at loading. The colors inside each UNF are proportional to their fissile isotope contents: red for high concentrations, blue for low and white for average (the fission product contents are not presented, but are usually inversely proportional to the fissile content). White and black rods indicate the position of the empty control rod guides. The canister geometry is schematically represented by the red circle (inside: light brown, outside: blue). As observed, each UNF has a unique pattern of fissile isotope concentrations, due to their unique irradiation life (e.g., reactor cycles, positions in the cores). These particularities play a role in the calculation of k\(_{\mathrm{eff}}\), as well as the relative position of each UNF. For these canisters, the k\(_{\mathrm{eff}}\) values are 0.94803 for the top one and 0.95533 for the bottom one [27]. Such difference is non-negligible as in the first case k\(_{\mathrm{eff}}\) is below the limit of 0.95, whereas in the second case it is above. This example illustrates the importance of the position of the UNFs within the canister. Note that also the orientation of each UNFs might influence the k\(_{\mathrm{eff}}\) value of the canister.

Fig. 2
figure 2

Examples of two canisters loaded with four different UNF assemblies (same ones are used in both canisters). Colors are proportional to the concentration of fissile materials (see text for details). For the top canister, k\(_{\mathrm{eff}}=0.94803\); for the bottom one, k\(_{\mathrm{eff}}=0.95533\)

Based on the above high resolution characterization of the UNF isotope contents, the decay heat as function of time is thus also computed by SNF for each individual material region. The decay heat value for a loaded canister is then simply obtained from a linear sum of the DH of the individual UNFs. These values do not depend on the relative location of the UNFs inside a canister and are computed by the SNF code (along with the inventories) [27, 28] at any given time (in this study, calculations are performed at 30 cooling/storage years). The k\(_{\mathrm{eff}}\) is calculated based on the same isotopic contents using a Monte Carlo neutronic simulation code (see next section), thus providing a consistent approach between the DH and k\(_{\mathrm{eff}}\) calculations.

3.2.2 Monte Carlo simulations for k \(_{\mathrm{eff}}\)

As mentioned, the Monte Carlo neutronic transport code used to compute k\(_{\mathrm{eff}}\) is called SERPENT [29]. Specific information is needed in order to perform such simulations: the canister geometry, its material, the location of the four UNFs (relative locations, axial and radial rotations), as well as the content of the assemblies: cladding, spacers and most importantly the actinides and fission products (basically actinides can generate neutrons, and fission products are neutron absorbers); the transfer of the isotopic compositions from the SNF code to the SERPENT code is realized with the help of the COMPLINK code (for details see Ref. [30]). In the current models, the different contents of each rod of each UNF are considered. Additionally, as mentioned, each UNF is typically divided vertically in 30 to 50 segments, the content of each segment being different due to the axial power profile in a reactor core. The canister geometry includes realistic dimensions, based on the model used in Ref. [27]. As it is customary in this type of safety analysis studies, the canister is considered flooded with water (apart from the UNF locations).

The SERPENT calculation time depends on the number of neutrons considered (also called neutron histories) and directly define the statistical precision of the calculated k\(_{\mathrm{eff}}\) values. In the criticality safety analysis, it is typical to perform calculations with a precision in the order of 10-50 pcm (or 0.01 %), and in the present study, the number of neutrons is selected to obtain a statistical uncertainty of about 10-15 pcm. As mentioned, the SERPENT calculations are usually time consuming and a first approximation is considered in the following by taking into account only a two-dimensional model with the central assemblies and canister axial segment (reflective boundaries are assumed). This allows to divide the calculation time by almost a factor 20. The consequence of this approximation is an artificial increase in the k\(_{\mathrm{eff}}\) values, but does not change the general study. Considering the 212 UNFs mentioned, various loading configurations have been randomly selected and a total of 46,000 k\(_{\mathrm{eff}}\) values were calculated. Figure 3 displays the number of appearances of k\(_{\mathrm{eff}}\) values given by SERPENT. The spread of the k\(_{\mathrm{eff}}\) values is relatively large, from about 0.75 to almost 1.05. Certainly such distribution depends on the selection of UNFs and on their random loadings, but it allows here to train a specific ANN and to explore its performances. One can notice that close to 30 % of the samples are above the safety requirement of 0.95. A similar distribution for the decay heat is presented in Fig. 4. As mentioned, such values are obtained by linearly summing the individual four components (one for each UNF), independently of their location and rotation in the canister. As observed, nearly 3% of the cases are also above the design limit.

Fig. 3
figure 3

Number of appearances of the k\(_{\mathrm{eff}}\) values from the \(\simeq\) 46000 Monte Carlo simulation samples computed after 30 years of cooling time

Fig. 4
figure 4

Number of appearances of the decay heat (DH) values from the \(\simeq\) 46,000 Monte Carlo simulation samples computed after 30 years of cooling time

This suggests that for the specific UNF designs investigated here, the DH might no longer be a limiting design factor for longer storage times (e.g., 40 years or more), consequently reducing the number of loading optimization constraints.

4 Development of an artificial neural network

As mentioned, artificial neural networks have been successfully used in multiple domains for regression problems. In the present study, the computing system is a supervised learning algorithm as it requires a labeled data set.

4.1 Description of the ANN

In order to predict the desired outputs (namely k\(_{\mathrm{eff}}\) values for specific canister loadings), the ANN needs to be trained. To that purpose, more than 46,000 samples (46 746 to be exact) were obtained with random canister loadings. In order to have an estimation of the spread of results due to the selection of training/testing sets, the 10-fold method is applied: the samples are once randomly mixed, then 10 % are used for testing the ANN and the rest is used for training. The procedure is repeated through the whole sample set, such that each sample is being used to test the ANN only once. Such approach would then provide 10 different biases (between the SERPENT k\(_{\mathrm{eff}}\) and the ANN k\(_{\mathrm{eff}}\)), providing a certain amount of confidence in the results, while keeping a reasonable training time (a few days). It is made with the Matlab environment and the “Neural Network Fitting app” [31].

The ANN’s loss function is the mean squared error function. The ANN is composed of fifty hidden layers; the training function is “trainbr”, using the Levenberg-Marquardt optimization to update the weight and bias values, and uses the bayesian regularization in order to find a neural ”network that generalize well” [31]. Other parameters are given in Table 1. Such selection of hyperparameters was performed using a trial and error process, and parameters leading to the best results are used in this study.

Table 1 Selection of parameters used for the ANN in the present work

4.2 Description of the input and output

As mentioned, each UNF is composed of 225 rods, with only 205 of them filled with fissile materials and fission products, the 20 other rods being empty. Initially, the SNF code provides 205 isotopes for each rod. As observed, the number of initial features is higher than the number of samples available, and there is therefore a need to reduce the dimension of the input. It is in particular motivated by the fact that some of these isotopes have a low concentration, or a low reaction rate, and their impact on the k\(_{\mathrm{eff}}\) is small enough to neglect them. The initial number of isotopes was then reduced to 38 by removing all isotopes below a certain concentration threshold, such threshold being also chosen by the trial and error process on k\(_{\mathrm{eff}}\). The remaining 38 isotopes are the following: \(^{87}\)Rb, \(^{90}\)Sr, \(^{90}\)Y, \(^{93}\)Zr, \(^{93,94}\)Nb, \(^{107}\)Pd, \(^{113}\)Cd, \(^{126}\)Sn, \(^{125}\)Sb, \(^{129}\)I, \(^{137}\)Cs, \(^{148}\)Nd, \(^{147}\)Pm, \(^{147,151}\)Sm, \(^{154,155}\)Eu, \(^{230}\)Th, \(^{231}\)Pa, \(^{233,234,235,236,238}\)U, \(^{237}\)Np, \(^{238-242}\)Pu, \(^{241,243}\)Am, \(^{242-246}\)Cm. From the aspect of nuclear cross sections, the remaining isotopes are indeed significant, given the characteristics of the UNFs.

Considering a limited number of 38 isotopes per rod, the total number of inputs for one specific k\(_{\mathrm{eff}}\) value is \(38\times 4\times 205=31,160\) values (4 UNFs, and each of them with 205 rods). Additionally, some isotopic concentrations are strongly correlated as concentrations are often slowly changing through radial neighbouring rods (see the smooth variations of colors in Fig. 2). It is therefore possible to reduce the input dimension, and a principal component analysis (PCA) is applied to the training set (and subsequently for all the 10 training sets from the 10-fold method). In the following, only the variables responsible for 99% of the variance will be used for the ANN, reducing the number of features to only 38. It was found that only 3300 features (inputs) are linearly independent from the 31,160 initial ones. Moreover, the case without PCA requires up to a few weeks to be computed and is thus not practical. Once performed, the PCA decomposition is stored and the same decomposition is applied to the testing set and with the genetic algorithm. All features were previously normalised with respect to the mean and standard deviation before the use of the PCA, otherwise only the highest concentrations will explain most of the variance. In the following, the “residual” is presented, defined as the difference between the k\(_{\mathrm{eff}}\) value computed by SERPENT and k\(_{\mathrm{eff}}\) value calculated by the ANN (expressed in pcm, or per cent mille).

4.3 Results of the ANN

After training the ANN with the 10-fold method, the mean residual for the test sets is 0 pcm, with an uncertainty of 1 pcm. Likewise, the standard deviation on the residual of the test sets is 45 ± 3 pcm. The quoted uncertainties are as one standard deviation (later mentioned as “std” or \(\sigma\)) given by the 10 ANN due to the 10-fold cross-validation.

As an example, Fig. 5 displays a comparison between the k\(_{\mathrm{eff}}\) values from SERPENT and from the ANN, for one specific case of the 10 test groups. This relation can be approximated by a linear regression. The R\(^2\) value of the training set is 0.999954, whereas it is 0.999932 for the testing set.

Fig. 5
figure 5

Comparison between the target value (k\(_{\mathrm{eff}}\) value from SERPENT) and the output (k\(_{\mathrm{eff}}\) from the ANN) for the testing set

In the following, only one ANN will be used with the genetic algorithm, being the one with the smallest standard deviation on the residuals for the test set (42 pcm). The bias of residual is 0 pcm and the extreme values are −740 and +467 pcm, see Fig. 6 for the distribution of the residuals.

Fig. 6
figure 6

Distribution of the residuals for the ANN, later used as an input for the genetic algorithm (1\(\sigma\) corresponds to 42 pcm)

As observed, the residual distribution is strongly peaked (with a kurtosis of 37.5), and only 1.6% of the cases are above 100 pcm (almost 2\(\sigma\)). It indicates that for the vast majority of the cases, the k\(_{\mathrm{eff}}\) values are extremely well reproduced by the selected ANN.

In order to find the best configuration, other designs for the neural network were tested. In Fig. 7, one can observe that results for neural networks with one or two layers. Also, two different activation functions for the hidden layers were investigated. It was then observed that the neural network with two layers and 50 neurons per layer using the sigmoid activation function provides the smallest standard deviation on the residual of the test set (45 ± 3 pcm). Additional configurations using three layers were also tested and did not show significant improvements. For clarity, these cases are not presented in the figure.

Fig. 7
figure 7

Standard deviation of the residual on the test sets depending on the number of neurons per layer, the number of layers and the activation function used in the hidden layer(s)

In order to assess the robustness of the model, the isotopic concentrations for each of the 38 isotopes of all the rods were perturbed using a normal distribution with a standard deviation of 1%. We found that the maximum difference between the nominal value from the ANN and the perturbed one was 40 pcm, corresponding to variations of the isotopic composition for \(^{238}\)U. Such value is still less than the standard deviation of the residual on the test set (42 pcm), as previously found.

4.4 Convergence

One should also ensure that the number of samples used to train and test the ANN is optimum. By choosing the same parameters as for the optimised ANN, the 10-fold procedure has been repeated with a reduced number of samples. The samples were still chosen randomly among the total set. Figure 8 presents the standard deviation on the testing set as a function of the total number of samples, used to both train and test the ANN. The last point is, therefore, the same value as found in Sect. 4.3. One can observe that the convergence seems to be reached for sample numbers above 40,000.

Fig. 8
figure 8

Standard deviation as a function of the number of samples. The last point corresponds to the one presented in Sect. 4.3

Such numbers certainly vary as a function of the initial UNF selection, and also as a function of how many UNFs are selected. But for the present study, the number of samples seems to be enough to reach convergence of the standard deviation.

5 Loading optimisation with the ANN and a genetic algorithm

Once the ANN is sufficiently trained, with a known residual and standard deviation, it can be confidently used in combination with a genetic algorithm (GA) to optimize the loading of canisters. Such final steps are presented in the following sections.

For the present problem definition, the GA is well suited as there is no specific structure for the “fitness function” to be optimized, but the search path is well constrained (minimisation of k\(_{\mathrm{eff}}\) and decay heat for each canister). Additionally, thanks to the ANN, the evaluation process is not computationally expensive.

5.1 Description of the genetic algorithm

The same previous UNFs are considered for the optimization of the canister loading: 212 UNFs to be loaded in a minimum of 53 canisters. As mentioned, their isotopic concentrations are determined from their irradiation history in the reactor core and cooling time, and are considered as fixed quantities in this study. For each canister, the first criteria are that their k\(_{\mathrm{eff}}\) and decay heat must be below the safety limits (0.95 for the criticality and 1.5 kW for the decay heat). The second criteria for the optimization are that the distribution of k\(_{\mathrm{eff}}\) and decay heat values for all canisters must be as homogeneous as possible. The GA will then pursue these four goals.

The genetic algorithm developed in that purpose is described in Fig. 9. At the initial step, all positions and orientations of each UNFs within the canisters are randomly chosen. The GA ends after 10,000 iterations. In each iteration, a loop is made over all canisters. In this loop, a random number will decide if the location or orientation of one UNF chosen will be randomly changed (within the same canister, or involving two canisters). If the k\(_{\mathrm{eff}}\) of the canister is smaller then the modification is accepted (if the swapping is within the same canister, the decay heat is not affected). In the case where the algorithm will swap an UNF with another canister, the alteration is kept if the maximum k\(_{\mathrm{eff}}\) and maximum decay heat between the two canisters are smaller than for the original configuration.

Fig. 9
figure 9

Schematic overview of the genetic algorithm for the canister loading optimisation. The k\(_{\mathrm{eff}}\) is evaluated with the trained ANN presented in the previous sections

5.2 Optimisation for k \(_{\mathrm{eff}}\) and the decay heat

Figures 10 and 11 display the evolution of the k\(_{\mathrm{eff}}\) and decay heat values as a function of the number of iterations in the GA.

As observed, the maximum k\(_{\mathrm{eff}}\) and decay heat values are well below the limits above 10 to 50 iterations. The final maxima are k\(^{\mathrm{max}}_{\mathrm{eff}}=0.9178\) and DH\(^{\mathrm{max}}=1.12\) kW (note that the k\(_{\mathrm{eff}}\) at iteration 0 is well above 1). Additionally, if the maximum values start to be unchanged above a few hundred iterations, the standard deviations are still decreasing up to a few thousand iterations (300 pcm for k\(_{\mathrm{eff}}\) and 11 W for the decay heat; the difference between maxima is 1800 pcm and 50 W for k\(_{\mathrm{eff}}\) and the decay heat, respectively). Taking into account the maximum absolute bias from the ANN (740 pcm, see Sect. 4.3) and the SERPENT statistical uncertainty (15 pcm) for k\(_{\mathrm{eff}}\), it can be confidently assumed that all safety criteria are fulfilled.

Fig. 10
figure 10

Maximum, minimum, mean and standard deviation of k\(_{\mathrm{eff}}\) for the 53 canisters as a function of the number of iterations while optimising at the same time k\(_{\mathrm{eff}}\) and the decay heat

Fig. 11
figure 11

Same as Fig. 10, but for the canister decay heat

5.3 Verification with SERPENT

Once the optimised loading pattern is found, one can check the agreement between the ANN-GA solutions and the SERPENT k\(_{\mathrm{eff}}\) values for the 53 canister loadings. To that purpose, additional 53 SERPENT calculations have been performed, corresponding to the loadings provided by the last iteration of the GA.

The agreement between the SERPENT k\(_{\mathrm{eff}}\) values and the ones from the ANN is found such that all the 53 k\(_{\mathrm{eff}}\) values agree within 740 pcm. Indeed, the maximum overestimation of the k\(_{\mathrm{eff}}\) by the ANN is of 340 pcm, and the maximum underestimation is 260 pcm. Apart from the cases at 340 and 260 pcm, all the other comparisons are within 120 pcm. The average deviation of the difference of k\(_{\mathrm{eff}}\) is 20 pcm. Given such small deviations, it can be concluded that the the ANN-GA results are validated by the reference SERPENT calculations.

6 Discussion of the present results

Based on these results, one can consider that the goal of our study is reached, namely optimizing the canister loading, by minimizing the k\(_{\mathrm{eff}}\) and decay heat values for each canister, as well as the the standard deviation of their distributions. Some additional comments are provided in this section.

6.1 Advantages of combining the ANN and GA

One of the goals of this study aimed at replacing the time-consuming Monte Carlo code SERPENT, the advantage being that once the ANN is trained, it only needs a fraction of a second to evaluate k\(_{\mathrm{eff}}\) for any specific canister loading (being an homogeneous or heterogeneous group of UNFs). Once the ANN trained, the use of a GA becomes possible. Such combination of ANN and GA allows to consider heterogeneous canister loadings from a practical aspect, opening additional optimization possibilities compared to existing approaches based on homogeneous loading (the direct use of SERPENT would inhibits the study of heterogeneous cases). For a more realistic case scenario with all the UNFs expected for the final repository, a new set of samples is required. Indeed, a large number of samples will be required to train sufficiently the ANN, which can also be limited by SERPENT’s computing time. The advantage of the ANN/GA combination is that the SERPENT calculation can be parallelized (random loadings), contrary to the SERPENT’s direct use in the GA.

6.2 Additional variability of the input data

The present study considers a number of 212 UNFs, with their fixed characteristics, coming from the “high resolution UNF simulations”. Such simulations are nevertheless based on a number of assumptions, such as nuclear data cross sections, manufacturing and engineering tolerances or impurities. These assumptions definitively impact k\(_{\mathrm{eff}}\) and DH values, as presented in Refs. [22, 28], and consequently can modify the training data for the ANN. It was for instance demonstrated in Ref. [33] that variation of nuclear data can lead to k\(_{\mathrm{eff}}\) distributions following the “Extreme Value Theory”, potentially heavily impacting the results presented in this study. Such variabilities are not considered in this work, as well as possible biases from the high resolution UNF simulations (for instance differences between the real isotopic concentrations and the calculated ones), and accountability errors (such as wrong UNF labelling or registration).

6.3 Current limitations and future improvements

The neural network is therefore a powerful tool to predict k\(_{\mathrm{eff}}\) values: if SERPENT calculations would have been run for 1000 GA iterations, more than 2000 days of Monte Carlo calculations would be needed (using a single CPU). However, its current implementation and the input data contain some limitations.

One limitation concerns the approximation of a two-dimensional geometry for the reference (SERPENT) calculations. Calculations with three-dimensional geometry can be performed, but it increases the number of input data (the multiplication factor is 40, as 40 vertical zones are considered). One needs to study the possibility of considering less zones, or apply a PCA with a limited impact on the k\(_{\mathrm{eff}}\) values computed by the ANN.

It is also expected that the current ANN is adequate for UNFs similar to the ones considered. Indeed, the ANN seems to reproduce the k\(_{\mathrm{eff}}\) values successfully for canister loaded with the 212 known UNFs. Nonetheless, the number of UNFs taken into account is small compared to the 12,000 used nuclear fuel expected for the final repository in Switzerland. The ability of the ANN to accurately compute k\(_{\mathrm{eff}}\) values for canisters with UNFs not yet considered during the training needs to be evaluated.

Moreover, under specific circumstances, one may wish to load a canister with fewer UNFs, if the decay heat of the UNFs is too high (for instance for the UNFs used in the last reactor cycle before decommissioning: such UNFs will have neutronic characteristics beyond the considered ranges). In that case, the current ANN needs to be re-evaluated.

An additional physical constraint not considered here is the availability of UNFs over time. Indeed, not all UNFs will be available at once at the docking/undocking stations of the final repository. A limited amount of transport casks can be unloaded at once, thus requiring a cask selection prior to the canister loading [21, 32]. Such practical limitation is not considered in this work. Additional aspects will need to be considered for the further development of the optimization methodology. Among other things, the following can be mentioned:

  • Dose rate constraints were not considered here while this constitutes a third key safety and design requirement. And when taking dose rates into account, the flexibility for the loading of the assemblies (e.g., rotation) might be reduced, meaning that the optimal loadings found here might no longer be valid. However, adding this third safety parameter will not change the results of the ANN to compute the k\(_{\mathrm{eff}}\) for the presented assemblies.

  • Only a static ANN/GA optimization was performed here since the k\(_{\mathrm{eff}}\) optimization was made at a given precise time, namely start of disposal. As mentioned earlier, the k\(_{\mathrm{eff}}\) might however vary with time and incorporating such dynamical effects in the ANN/GA optimization scheme will need to be assessed.

  • The usage of high resolution UNF simulations combined with the assumption that UNFs could be loaded in any arbitrary manner (e.g., orientation, rotation) was shown to allow using the exact minimum number of required canisters. However, it remains to be verified if this could also be achieved (or if a higher number of canisters would be required) when the optimization is performed without allowing anymore for arbitrary assembly placements and/or when assuming assembly-average UNF contents. Addressing this is relevant since as shown here, arbitrary placements could be advantageous from a loading point of view. On the other hand, it could reduce flexibility in terms of manageability and operations. And it would also lead to stricter requirements in terms of safety assessments since the risk for inadvertent fuel assembly loadings due to machine and/or human errors would become larger.

In the light of these limitations, additional efforts are needed towards the full application of this method for industrial goals. A large part of the necessary information is already in place (realistic input data), but its full usage for the optimization of the canister loading is not yet fully demonstrated.

7 Conclusion

This work is oriented towards finding optimised canister loadings for used nuclear fuel for the Swiss deep geological repository. The optimization was performed on two quantities: the canister criticality and its UNF decay heat. Additionally, uniform distributions for both quantities were required. This study contains a description of the surrogate model used to replace the Monte Carlo transport code to compute the k\(_{\mathrm{eff}}\) for each canister and presents the genetic algorithm applied to optimise the loading.

The surrogate model is an artificial neural network aiming at evaluating k\(_{\mathrm{eff}}\) values of canisters. The artificial neural network was therefore trained and tested on a limited set of UNFs. The evaluation of the surrogate model on the testing set used was presented. The test set was composed of the same UNFs as in the training set but loaded differently.

Once a surrogate model has been developed, a genetic algorithm was used to optimise the loading, based on the ANN evaluations. The same UNFs used to train and test the model have been used for this loading optimisation problem. The final optimized loading complied with the requirements, namely that all quantities were below maximum allowed values, and that the distributions of decay heat and k\(_{\mathrm{eff}}\) among canisters were as uniform as possible. Finally, a number of points of improvements were presented, but the good performances of the present results are encouraging for future developments, involving more realistic constrains.