This study considers a 17 × 17 PWR fuel assembly as the reference design, consisting of 264 fuel rods of Fuel Type 1 with UO₂, 32 fuel pins of Fuel Type 2, consisting of UO₂ mixed with burnable absorber, and 25 guide tubes filled with water, as shown in Fig. 2. Table 1 lists the detailed specifications of the ref. design. This study was conducted under cold zero power (CZP) conditions, maintaining constant temperatures of 293K for fuel, cladding, and gas gap.

A total of 06 (six) decision variables were considered for training the ANN model:

The reason for selecting these decision variables is that they all influence the value of the output parameter, i.e., infinite multiplication factor (k_∞) of an assembly, as shown in Fig. 3. These decision variables were randomly varied in the range specified in Table 2 to obtain corresponding values of k_∞ from OpenMC. Uniform distribution was considered for the random datapoints to avoid sample bias in the dataset. Thus, the trained ANN model is continuous in nature, meaning it can give continuous values of k_∞, unlike the other ML models in the literature that give discrete outputs like classification problems (Akkurt and Çolak 2002; Cai et al. 2021). The materials and their corresponding densities used in the OpenMC simulation are outlined in Table 3.

The range of U-235 enrichment and burnable poison contents was selected so that the ANN models are trained for BOC to EOC conditions of the assembly. Additionally, the water density varied such that the ANN can simulate accident scenarios with very high void fraction. Finally, the concentration of soluble poison in the coolant was also taken as a variable, keeping in mind the need for negative reactivity insertion in case of an emergency.

Data-driven models require extensive data to effectively learn and analyze patterns to establish relationships between input parameters and the target output. In this study, data points were simulated using OpenMC version 0.13 and the ENDF/B-VIII.0 cross-section library. OpenMC employs random sampling techniques and the stochastic movement of neutrons within a stochastic process that mimics a real experimental framework. The process flow of Monte Carlo Neutronic simulation is illustrated in Fig. 4. This study did not focus on burn-up analysis, thus only a cross-section data library was required and not the depletion chains.

The Monte Carlo simulations were performed for a steady-state reactor operating condition. Furthermore, reflective boundary conditions were considered in this study to get the value of infinite multiplication factor. A total of 6000 OpenMC datapoints were generated from which 70%, 15% and 15% were utilized for training, testing, and validation, respectively. Each simulation was executed with 200 batches, including 20 initially inactive batches, each batch consisting of 40,000 particles per simulation. These simulation parameters resulted in an uncertainty of less than 50 pcm.

An Artificial Neural Network is an ML model that mimics the architecture and functionality of the human brain. It learns by adjusting synaptic connections, often referred to as updating weights and biases. It consists of multiple layers, which are interconnected with nodes, where each performs a basic calculation to produce an output. Next, the output is passed sequentially from one layer to another, and the final output of the architecture serves as either a prediction or classification. Typically, ANNs consist of three main components: input layer, hidden layer, and output layer. All training data are supplied to the input layer, processed in the hidden layer, before reaching the output layer. The complex architecture of ANN makes it a powerful tool that has been successfully utilized in numerous research fields (Abiodun et al. 2018; Ghritlahre and Prasad 2018; Ghani et al. 2022). Fig. 5 shows a typical four-layer ANN architecture. Each layer in the architecture consists of interconnected nodes, or neurons, each performing basic computation on its inputs. The input to a neuron is multiplied by a specific weight and then summed together. The resulting sum is then passed through an activation function, introducing non-linearity to the neuron’s output. This output becomes the input for the neurons in the next layer.

The weights and biases within the ANN architecture are optimized using the backpropagation technique. During training, the network is fed with input data labelled with the true output. Then the weights and biases are updated to minimize the discrepancy between the true output and the predicted output. The training process of ANN model is illustrated in Fig. 6.

The output in a node and the net output are calculated as,

yout =b+wixiT (1)

ynet =E(yout ) (2)

Here, x, w, and b represent the inputs, weights, and bias. An activation function ‘↋’ is imposed on the node output and introduces non-linearity. The x, w is expressed as follows.

x=[ x 1 x 2 x 3⋮ x n] and w=[ w 1 w 2 w 3⋮ w n] (3)

The number of input parameters is represented by X and the characteristics of the output node Y_net are determined by the non-linear activation function ↋(.).

Some of the most common activation functions which are frequently used in ANN are summarized in Table 4. The choice of activation function depends on the area to be applied. In this study, the optimal activation function for the hidden layers was evaluated by comparing the training mean squared errors (MSE). The Levenberg-Marquardt algorithm was employed during the training. This work employs the Neural Network Fitting Toolbox in MATLAB 2018 for developing the ANN models.

This study considered ANN models consisting of one input layer with six nodes for the decision variables, one or two hidden layers with a varied number of nodes, and an output layer with two nodes. The reason behind two output nodes is the inherent uncertainty of Monte Carlo simulation, which results in an upper and a lower limit of k.

Two distinct ANN architectures were considered:

The authors also trained two 5 layer ANNs but as a supplementary check. The maximum number of neurons in each hidden layer was set based on the rule of thumb that states the training data should be 10 times the number of parameters of the ANN. Since the available dataset consisted of 6,000 points, the upper limit of Z, X and Y were set to 50, 20 and 15, respectively.

While the rule of thumb is to supply training data ten times the number of neurons, a recent research recommends increasing this factor to 50 times, which can yield superbly accurate ANN models (Alwosheel et al. 2018). In line with this, the present study rigorously optimizes the ANN models by varying the number of neurons in the hidden layer. The optimization process was initiated by evaluating the most suitable activation function. Firstly, a 3-layer ANN architecture with 20 hidden neurons and a 4-layer architecture with 10 neurons per hidden layer were considered. Next, the number of hidden layer neurons was increased and decreased to determine the optimal configuration, while the previously selected activation function was kept unchanged. During the training process, a maximum number of epochs was set to 1000 to prevent the possibility of overfitting of the ANN models. Additionally, the training was stopped if the validation MSE remained unchanged for six consecutive epochs.

All ANN models were trained in MATLAB 2018a software with a dataset of 6,000 datapoints obtained from the OpenMC simulations. The dataset was divided into three subsets for training (70%), validation (15%) and testing (15%), respectively. The training performance of each model was evaluated based on its training MSEs given by,

MSE=∑i=1N(yi−y^i)2N (4)

Where y_i and ŷ_i represent actual and predicted output results, respectively, and N is the total number of data points. The models that satisfied the acceptance criteria of training MSE ≤ 25 × 10^-6 were further tested with a separate dataset of 117 points to ensure they have predictive errors of ≤ 500 pcm.

Void fraction (α) and void coefficient (α_v) are important parameters in reactor core design and safety management. The objective of this step is to investigate the impact of void fraction on the void coefficient and predict reactivity changes to ensure the reactor’s normal operating conditions at required material compositions. These aspects are essential for the design optimization of a 17 × 17 PWR fuel assemblies. The expression for void fraction is,

α=11−ρgρf(1−xx) (10)

Here ρ_g and ρ_f are densities of saturated steam and saturated liquid, respectively, and x denotes the quality of coolant defined as,

x=(v−vf)vfg (11)

Where values found from the steam table at 1 atm are: the specific volume of saturated vapor, ν_g, is 1.6734 m³/kg, and the specific volume of saturated liquid, ν_f, is 0.001043 m³/kg. And the expression void coefficient is,

αv=δρδα (12)

Where the reactivity is defined as,

ρ=(k−1k) (13)

Reactivity is sensitive to water density; depending on the moderator-to-fuel ratio, reactivity may either increase or decrease with a decrease in water density. In the reactor, accidents leading to boiling crisis and bubbles (or voids) formation at higher temperatures will decrease water density. Therefore, it is desired that the reactor power, which is a function of neutron flux, goes down as the void volume increases. This is possible if the reactor core has a negative void coefficient, which means the reactivity will decrease with an increasing void fraction. A negative void coefficient is always preferred to control reactivity, prevent super-criticality, and maintain normal operating conditions of the reactor core. It also enhances safety by serving as a built-in safety mechanism in a water-moderated reactor core, improving reactor controllability under adverse conditions such as coolant loss (Ishida et al. 2023; Jandaghian et al. 2024).

A negative void coefficient is achieved if the assembly is under-moderated, which means it is to the left of the moderator-to-fuel ratio that corresponds to the maximum value of k_∞ of that assembly. If the assembly pitch is selected properly, it is possible to make sure that the void coefficient remains negative throughout the operating lifetime of a fuel assembly. Thus, the main purpose of this step is to optimize the 17 × 17 PWR design by finding a fuel pitch that ensures the reactor’s under-moderated status from BOC to EOC conditions.

Abbreviations:

ANN Artificial Neural Network;

BOC Beginning-of-cycle;

CZP Cold Zero Power;

EOC End-of-cycle;

GA Genetic Algorithm;

GPR Gaussian Process Regression;

MAE Mean Absolute Error;

ML Machine Learning;

MSE Mean Square Error;

pcm Percent milli-rho;

PSO Particle Swarm Optimization;

PWR Pressurized Water Reactor;

SA Simulated Annealing;

SS Simplex Search;

SVR Support Vector Regression;

TANH Hyperbolic Tangent.

Symbols:

α Void fraction;

α_v Void coefficient;

ϵ Fast Fission Factor;

η Reproduction Factor;

ν_f Specific volume of saturated liquid (m³/kg);

ν_g Specific volume of saturated vapor (m³/kg);

ρ Reactivity (pcm);

ρ_f Saturated liquid densities;

ρ_g Saturated steam densities;

f Thermal Utilization Factor;

k_∞ Infinite Multiplication factor;

N_m/N_u Moderator-to-fuel ratio;

OF Objective Function;

P Resonance Escape Probability;

P_t Fuel Pitch.

This study considers various activation functions for all possible combinations of two ANN architectures: 6-Z-2, 6-X-Y-2. Results revealed the tanh and sigmoid activation functions exhibit superior performance for 3-layer and 4-layer ANN models, respectively (Suppl. material 1: table SS1). In contrast, ANN models which used ReLU in the hidden layers resulted in the highest training MSEs, indicating poor performance. As shown in Fig. 7(a), increasing the number of hidden neurons beyond 40 did not yield any further improvements in MSE. Additionally, 3-layer ANN architectures (except 30, 40, and 50 hidden neurons) failed to meet the acceptance threshold of 25 × 10^-6. In contrast, several 4-layer ANN architectures successfully satisfied the acceptance criteria, with the 6-15-15-2 model achieving the lowest testing MSE of 0.98 × 10^-6. The MSE of this architecture was also lower than 3-layer architecture (Fig. 8(b)).

Although 6-X-Y-2 ANN models with X ≥ 10 satisfied the acceptance criteria, when the training and testing MSEs are compared, it is evident that the 6-15-15-2 ANN model has the lowest training MSE among all other architectures. This is further illustrated in Fig. 8, which clearly shows the minimal difference between training and testing MSEs for the 6-15-15-2 architecture (Fig. 8(a)), while a significant discrepancy is observed for the 6-15-10-2 ANN model (Fig. 8(b)). Both models demonstrated a steady decline in training, validation, and testing MSEs up to approximately 400 epochs, which indicates a slight underfitting scenario rather than overfitting. The curves exhibited a flattening trend, suggesting that the ANNs had nearly converged to their optimal performance. Since the final MSE values were below the threshold, extending the training further was redundant. Although two 5-layer ANNs were also trained, their performance was unsatisfactory, with greater deviation between training and testing MSE (Suppl. material 1: table SS1). Based on these findings, the study finalized this 4-layer ANN architecture with 15 and 15 neurons in the hidden layers along with sigmoid activation functions for efficiently predicting neutronic parameters of a 17 × 17 PWR fuel assembly.

While ANNs are known for their capabilities in learning very complex correlations, some studies have shown the superior prediction performance of Support Vector Regression (SVR) and Gaussian Process Regression (GPR) over the ANNs and vice versa (Jia et al. 2019; Lawal et al. 2023), which is why this work considered them and a comparative analysis was performed (Suppl. material 1: table SS2). It is obvious from the table that the ANN model outperformed both GPR and SVR models, where neither GPR nor SVR models were able to achieve the target MSE of 25 × 10^-6. The optimized ANN architecture (6-15-15-2) achieved the lowest testing MSE of 0.98 × 10^-6. In contrast, GPR models, although showing comparatively better performance than SVR models, still resulted in errors two to three times the ANN’s MSE. The Quadratic kernel yielded an MSE of 135.4 × 10^-6, best among all variants. The SVR models performed the worst, with the lowest error found for Cubic kernel with an MSE of 661 × 10^-6, which was beyond the threshold of 25 × 10^-6. Based on the results, it is opined that the ANN architecture (6-15-15-2) is more capable of handling the complex dynamic relationship among design variables for a 17 × 17 PWR fuel assembly compared to GPR and SVR. Therefore, ANN is one of the most reliable and efficient frameworks for developing an efficient model in this context.

The regression plots for the developed ANN architecture are illustrated in Fig. 9. It shows the goodness of fit for the ANN model. Here, only 1500 random data points were shown for better visualization. As noted earlier, the dataset was divided into 4200 training points, 900 validation points, and 900 testing points. The overall R² values for this data-driven model were found to be close to unity. It indicates an excellent agreement between predictions and reference data. Thus, it is obvious that the developed ANN model exhibits superb accuracy.

Fig. 10 illustrates the error histogram for the developed ANN model. It is obvious that the error distribution exhibits a near-Gaussian (bell-shaped) pattern, where a significant portion of errors is condensed around zero. The spread of errors beyond ± 3.0 × 10^-3 was relatively rare, indicating that large deviations were minimal. Overall, the predictive error was symmetrically distributed and remained well within the acceptable accuracy threshold, confirming the robustness of the ANN model. It is noteworthy that only 91.3 pcmMAE was obtained from the error histogram (Fig. 10) for the developed ANN architecture, which is below the threshold of 500 pcm.

Table 6 shows the results obtained from the sensitivity analysis of the developed ANN model. From Table 6, it is clear that all the input parameters fulfill Condition 3 in Table 5, indicating that they have non-linear correlation with the output parameters. However, k_∞ seemed to be comparatively less sensitive to UO₂ density, which was somewhat expected considering that the small variation of fuel density in the training set (0.97 g/cm³) shouldn’t have caused a large change in multiplication factor. On the other hand, the output parameters were found to be highly sensitive to the concentration of burnable poison (Gd₂O₃) in the fuel and moderator density, both agreeing with the underlying physics. Finally, k_∞ was sufficiently sensitive to the other three input parameters, indicating that the choice of the decision variables for the ANN model was appropriate. It is noteworthy that the sensitivity values (S¯ik±Sikσ) were almost identical for k_{∞, max} and k_{∞, min}, which was completely expected since the two output values are nothing but the upper and lower limits the same variable, i.e., multiplication factor.

The performance of the developed 6-15-15-2 ANN model confirmed that it can efficiently mimic the complex, non-linear relationship between the input and output variables for predicting the neutronic parameters of 17 × 17 PWR fuel assemblies. Although it is a pure mathematical model without any physics-based constraints, it showed an inherent ability to capture the intricate relationship among various variables of 17 × 17 PWR fuel assemblies and neutronic parameters. Therefore, all observations suggest that the developed ANN model is extremely reliable for the precise prediction of the neutronic parameters of PWR fuel assemblies. Moreover, findings from the sensitivity analysis suggest that the ANN model is well-fitted, as an overfit model would have revealed unphysical characteristics.

The prediction accuracy of the ANN model was re-tested with a separate testing dataset generated from OpenMC, since Monte Carlo simulations are considered the “gold standard” for neutronic analysis (Horelik et al. 2014). Fig. 11 presents a comparison between ANN-predicted and OpenMC-obtained values of k_∞ for the reference geometry of 17 × 17 PWR assembly. Remarkable agreement between OpenMC and ANN models is observed on Fig. 11. The average predictive errors for the ANN model were around 73.77, 71.98, and 105.73 pcm for 1, 3, and 5 wt% Gd₂O₃ content, respectively. These values are way below the threshold of 500 pcm, establishing that the ANN models are sufficiently accurate.

From Fig. 11, the typical relationship between k_∞ and P_t is also observed, where the multiplication factor initially increases with fuel pitch in the under-moderated region but starts decreasing as the assembly becomes over-moderated. For visualization purposes, the curves join the maximum k_∞ values for different U-235 enrichment levels. It may be observed that the curves resemble straight lines, meaning that the selection of multivariate regression for correlation development is appropriate.

To enhance the computational efficiency of the developed 6-15-15-2 ANN architecture, four different optimization algorithms: Simulated Annealing, Simplex Search, Particle Swarm Optimization, and Genetic Algorithm were tested under nine distinct simulation cases. These cases accounted for variations in U-235 enrichment levels (3, 4, and 5 wt%) and Gd₂O₃ contents (1, 3, and 5 wt%), while the values for UO₂, B10 and Water were kept constant at 10.41 g/cm³, 0 ppm and 0.74 g/cm³.

The results demonstrate that SS-based simulation coupled with the developed 6-15-15-2 ANN architecture resulted in a significant computational advantage across all tested algorithms. As shown in Table 7, the ANN architecture had faster convergence when coupled with SS and PSO, with average computation times of 0.646 s and 0.783 s, respectively. In comparison, GA required significantly longer time of 2.158 s, while SA performed the worst, taking 4.189 s, which is almost twice as much as GA, highlighting its inefficiencies. These findings indicate simplex search is the most efficient optimization tool under the given constraints.

It is important to note that SA required a higher computational time compared to GA and PSO. This was primarily due to SA’s inability to meet the target tolerance even after 800 iterations, along with its continued improvements beyond 50 consecutive iterations. However, when the time per iteration is considered, SA was faster than both GA and PSO.

The simplex search algorithm, coupled with the developed ANN architecture, was utilized to determine P_t that corresponds to k_{∞, peak} for different U-235 enrichment levels and Gd₂O₃ concentrations. The data obtained from SS simulations were then utilized to derive suitable correlations for predicting P_t and k_{∞, peak} using multivariate linear regression. The equations derived are presented in Table 8, both having R²-values greater than 0.90. From the equations, it is clear that k_{∞, peak} increases with increasing U-235 enrichment (w_U235) while decreasing with Gd₂O₃ content (w_Gd2O3). Assuming EOC U-235 enrichment as 1.0 wt% and complete burning up of the burnable poisons, this study suggests that the fuel pitch should be ≤ 1.434 cm to keep the assembly under-moderated throughout its operational lifetime, ensuring neutronic safety. For a BOC assembly with maximum possible U-235 enrichment (5.0 wt%) and 3.0 wt% Gd₂O₃ content, the assembly will be under-moderated if fuel pitch is ≤ 1.6749 cm, which shows the importance of considering the neutronic behavior of both BOC and EOC fuel assemblies for proper evaluation of safety margin.

Figs 12, 13, 14 show the changes in k_∞ and α_v with α for three different fuel pitches: 1.5 cm, 1.8 cm, and 2 cm, respectively, all under CZP conditions and in the absence of soluble poison in the coolant. The fuel enrichment, burnable poison concentration and water density were varied while UO₂ density (10.41 g/cm³) and boron concentration (0 ppm) were kept constant. From Fig. 12, it is observed that the void coefficient of the assembly with 1.5 cm fuel pitch may be positive at the initial stages of void formation, depending on U-235 enrichment and Gd₂O₃ concentration. Increasing the concentration of burnable poison and decreasing enrichment is needed to achieve a negative void coefficient for all possible void fractions. From Figs 13, 14, however, it is evident that the higher fuel pitch will result in positive void coefficients even at α ≤ 0.35, which is completely unsafe for the reactor. Therefore, the fuel pitch should be kept smaller than 1.5 cm to ensure the inherent safety of the fuel assembly and thus, the reactor core. This observation is consistent with the results from the derived equations (P_t ≤ 1.434 cm).