Introduction

With the depth of its intellectual development, mathematical modeling and simulation (M&S) have emerged as a powerful tool that promises to revolutionize engineering and science. M&S are expected to describe a system, on the basis of which further research can be conducted. For example, the Consortium for the Advanced Simulation of Light Water Reactors (CASL), was established by the US Department of Energy (DOE) in 2010, with the goal of providing M&S capabilities that support and accelerate the improvement of nuclear energy’s economic competitiveness and ensure nuclear safety [1].

Numerical simulations can differ significantly from experimental observations, and minimizing this difference has always been an important topic in engineering practice [2-4]. The difference between model and reality is particularly true in reactor physics, primarily for the following reasons:

(1) The complex nature of nuclear phenomena. A nuclear system may involve neutron transport, thermal hydraulics, and fuel performance, among others, which makes it difficult to accurately model the neutronic behavior. Theoretically, this can be well described by the Boltzmann transport equation, incorporating coefficients derived from various experimentally evaluated neutron interaction cross sections. These cross sections have strong, discontinuous behavior in space due to material heterogeneities and extreme variations with energy due to resonance phenomena associated with compound nucleus formation [5-7].

(2) Inaccuracy of input parameters. A model requires input design data or explanatory data, such as geometry, composition, and control parameters, and input model parameters, including physical constants such as interaction cross sections and material thermal and mechanical properties. Particularly, the two-level approach (i.e., a micro-level model for determining macro cross sections and a macro-level model for determining the final response of interests) for reactor physics calculation may lead inexact input parameters [5-7]. As a consequence, inexact parameters will clearly cause inaccurate outputs.

With the advancement of science and computational power, we can probably make fewer simplifications and approximations in our models, resulting to reduced differences, whereas the resources spent must still be considered for economic reasons. As we reduce this difference, it becomes increasingly dominated by the input parameters.

In particular, nuclear data, such as neutron cross sections and resonance parameters, are evaluated by combining approximated nuclear physics and experimental observations, which have epistemic and aleatoric uncertainties. These uncertainties then propagate through a neutronics model, creating differences and uncertainties [8]. It is easy to imagine that it would be the best option if we were able to re-evaluate the nuclear data. Unfortunately, the nuclear evaluation process is significantly long and costly, making it unrealistic. For example, [9] estimated the cost of a single new observation of a nuclear datum at 400,000 US dollars.

The uncertainty of the input parameters causes uncertain behavior in the model, resulting to the development of three fields in statistics: uncertainty quantification (UQ), sensitivity analysis (SA), and data assimilation (DA) [10]. They are beneficial in engineering design, with the potential to significantly reduce unnecessary costs. For example, knowledge of the sensitivity of attributes has been used to reduce costly design margins.

Data assimilation assimilates previous experience or experimental observations to reduce the differences and uncertainties in model predictions. It was first proposed in the field of meteorology and is widely used to predict meteorological phenomena, such as hurricanes [11]. Data assimilation methods can be roughly divided into two categories: variational data assimilation [12-14] and sequential data assimilation [15-17], whereas some hybrid data assimilation techniques have emerged [18-20]. Bayesian approaches to data assimilation were established after [21, 22] and [23] provided an overview of the Bayesian perspective, discussing some approaches from this perspective. In neutronics, data assimilation is primarily used to predict the bias of critical systems [24], thereby reducing the uncertainty of advanced reactor designs [25]. In recent years, research interest in data assimilation in neutrons has mainly focused on the adjustment of nuclear data, which could establish cross-correlations between nuclear data or nuclides that are traditionally unavailable [26, 27]. In addition, new methods based on stochastic sampling of input parameters have been proposed, such as MOCABA [28] and BMC [29]. In addition, D. Siefman showed that these two methods yield similar results to each other, as well as the traditional method known as generalized linear least squares (GLLS) [30]. Readers can refer to [31-33] for the Bayesian applications in thermal hydraulics.

The method proposed in this study involves sampling the probability distribution, and the method we used is MCMC. The MCMC method began with the Metropolis method proposed by Metropolis et al. in 1953 and was extended byW.K. Hastings in 1970. In 1984, S. Geman et al. demonstrated how the method, known as the Gibbs sampler, can be adapted to high-dimensional problems that arise in Bayesian statistics [34]. In 1987, S. Duane et al. combined MCMC and molecular dynamics methods, known as the hybrid Monte Carlo method or Hamiltonian Monte Carlo (HMC) [35]. Later, scholars developed a series of improved HMC, such as the No-U-Turn Sampler (NUTS) [36] and stochastic gradient Langevin dynamics (SGLD) [37], based on HMC.

Although the Bayesian approach has been widely used in data assimilation, there is still a lack of research on the selection of a prior distribution, which is generally selected as the probability representation that exists before experimental observations, and is probably unsatisfactory. It is evident that the prior distribution affects the results of the method. This study aims to complete this task, namely, to provide the initial input of prior distributions for classical data assimilation methods. As an example, recent work [38] has applied the Bayesian neural network method to predict the beta decay lifetime of atomic nuclei and their uncertainty, which improves learning accuracy and uses prior distributions of parameters; therefore, it may be useful to consider using the method proposed in this work once before the algorithm to improve the performance of the method with more accurate prior distributions. Moreover, most data assimilation methods only provide point estimations for parameters; however, what if nuclear data are actually probabilistic? In fact, certain parameters in our understanding resemble distributions or random variables, such as the number of neutrons released per fission and the change in the concentration of nuclear fuel over a short period due to fuel depletion. We attempt to fill this research gap and take the first step in this distribution inference. Therefore, we propose an ensemble Bayesian method. It is also worth mentioning the following:

(1) The proposed method is simulation-based. As strategies for the design and evaluation of complex nuclear systems have shifted from heavy reliance on expensive experimental validation to highly accurate numerical simulations, more attention has been paid to simulation-based approaches. This allows our method to play a role in guiding the overall system design as well as optimizing the setup of the validating benchmark experiments.

(2) The proposed method has good compatibility. Because there are numerous legacy codes used extensively in the design process, they are expected to persist as integral parts of nuclear design calculations. Our method can be incorporated in the legacy codes, making it tractable for nuclear systems simulation.

The remainder of this paper is organized as follows. In Sect. 2, we describe the parameter distribution inference problem in a general manner and describe the corresponding experiment. In Sect. 3, we first introduce the Bayesian inference and Metropolis-Hastings algorithms, followed by the core algorithm. In Sect. 4, we analyze the convergence of the method and illustrate a series of numerical tests on the proposed method, including a finite cylindrical reactor and the 2D IAEA benchmark problem. Finally, we provide a brief conclusion in the last section.

Problem Description

When modeling a physical phenomenon, the first step is to introduce reliable theories, such as conservation laws and transport theory. In nuclear reactor physics, there are many classic models such as the finite cylindrical reactor and the 2D IAEA benchmark problem, which will also be presented in Sect. 4. However, it is worth emphasizing that the method proposed in this study is valid, regardless of the theories introduced. For clarity, we skip this step and consider the following general system: $\{\begin{array}{l}y(x\mathrm{,}\mu )=f(x\mathrm{,}\mu )\hfill \\ z(x\mathrm{,}\mu )=y(x\mathrm{,}\mu )(1+ϵ)\hfill \end{array},$ (1) where $x\in D\subset R^3$ is a spacial vector, $y\in R$ is an observable physical quantity, $z\in R$ is an observation of y possibly with noise ϵ, and $\mu ={\left({\mu }^{\left(1\right)}\mathrm{,}\cdots \mathrm{,}{\mu }^{\left(n\right)}\right)}^{\text{T}}$ are input parameters of the model.

f represents the intrinsic functional relationship between these variables, which models the corresponding physical phenomenon. $ϵ\sim N\left(\mathrm{0,}{\sigma }^2\right)$ is the zero-mean Gaussian relative noise that models the observation noise.

Ensemble Bayesian Method

The Bayesian approach has been widely used in data assimilation and is the method proposed in this study. Subsequently, we provide an overview of the Bayesian inference.

Numerical Examples

In this section, we examine the effectiveness of the proposed method. Two models in nuclear reactor physics, namely the finite cylindrical reactor and the 2D IAEA benchmark problem, are tested below. The first test case, adapted from [42], has an analytical solution; the 2D IAEA benchmark problem is one of a classical benchmark problem in nuclear physics [42], and has only a numerical solution. To verify the effectiveness of the distribution estimation method, some parameters in these models were assumed to be uncertain in the following tests:

4.1

The Finite Cylindrical Reactor

A finite cylindrical reactor is a multiplying system of a uniform reactor in the shape of a cylinder of physical radius R and height H. Its behavior is modeled using the following two-dimensional monoenergetic diffusion equation: $\{\begin{array}{l}D\Delta \varphi \left(r\mathrm{,}z\right)-{\text{Σ}}_{\text{a}}\varphi \left(r\mathrm{,}z\right)=\frac{1}{k_{\text{eff}}}\nu {\displaystyle \sum {}_{\text{f}}}\varphi \left(r\mathrm{,}z\right)\hfill \\ \varphi \left(R_{\text{e}}\mathrm{,}z\right)=\mathrm{0,}\varphi \left(r\mathrm{,}{H}_{\text{e}}\right)=\mathrm{0,}r\in \left(\mathrm{0,}{R}_{\text{e}}\right)\mathrm{,}z\in \left(\mathrm{0,}{H}_{\text{e}}\right)\hfill \end{array},$ (16) where $R_{\text{e}}=R+d\mathrm{,}{H}_{\text{e}}=H+d$, and d is the extrapolated length.

By replacing the Laplacian with its cylindrical form, as shown in Fig. 1, the analytical solution of (16) is as follows when k_eff=1: $\varphi \left(r\right)=A\cdot J_0\left(\frac{2.405r}{R_{\text{e}}}\right)\cos \left(\frac{\pi z}{H_{\text{e}}}\right)\mathrm{,}$ (17) where J₀(r) is a zeroth-order first-order Bessel function.

Fig. 1Full-screen View

(Color online) Geometry of the finite cylindrical reactor [43].

Our goal is to estimate the distribution of certain uncertain parameters from the observations. Later, we compare two alternative schemes and examine the convergence of the method by applying them to single-parameter estimates, followed by multiparameter estimates. We set A=20,H_e=18 and use the ensemble Bayesian method proposed in the previous section to estimate the distribution of the parameter R_e. In other words, we consider the following system: $\{\begin{array}{l}y\left(x\mathrm{,}\mu \right)=20\cdot \left(\frac{2.405x}{\mu }\right)\cos \left(\frac{\pi x}{18}\right)\hfill \\ z\left(x\mathrm{,}\mu \right)=y\left(x\mathrm{,}\mu \right)\left(1+ϵ\right)\hfill \end{array},$ (18) where $x\in [\mathrm{0.5,7.5}]\subset R$ is one-dimensional spatial position, μ is the input parameter we are interested in and is assumed to have a real distribution, J₀ is the zeroth-order first-order Bessel function. $y\in R$ is an observable physical quantity, and $z\in R$ is the observation of y with noise ϵ. In this study, observations were prepared using synthetic data from simulations rather than real observations.

To prepare the observation data for our experiments, we first select m positions at equal intervals in [0.5,7.5] and repeatedly obtain the observations for the corresponding positions under k pairs of μ sampled from the real distribution. Suppose $ϵ\sim N\left(\mathrm{0,}\sqrt{0.2}\right)$, the real distribution of μ is N(11.0,0.36) and the prior distribution is N(12.0,0.25), which could be an empirical estimate. In this case, the number of sensors m, number of samples k, number of convergence steps, and total number of steps of the Metropolis-Hastings algorithms N_con and N_tot, which are introduced in Sect. 2.2, affect the results. Subsequently, we fix the value of $N=N_{\text{tot}}-N_{\text{con}}=100$ and discuss the effect of $m\mathrm{,}k\mathrm{,}{N}_{\text{con}}$ on the result. Before doing so, we compare two alternative schemes after sampling from the posterior distribution $p\left({\mu }_j|z_{\mathrm{1,}j}\mathrm{,}\cdots \mathrm{,}{z}_{m\mathrm{,}j}\right)$.

4.1.1

Comparison of Two Alternatives Schemes

In this subsection, we compare these two schemes. Because ${\mu }_1\mathrm{,}\cdots \mathrm{,}{\mu }_k$ are k samples of μ, a straightforward idea is to restore ${\mu }_1\mathrm{,...,}{\mu }_k$. After introducing the Metropolis-Hastings algorithm, we can sample from the posterior distribution $p\left({\mu }_j|z_{\mathrm{1,}j}\mathrm{,}\cdots s\mathrm{,}{z}_{m\mathrm{,}j}\right)$ and traditional practice is to estimate μi using samples from $p\left({\mu }_j|z_{\mathrm{1,}j}\mathrm{,}\cdots \mathrm{,}{z}_{m\mathrm{,}j}\right)$. Another scheme is to aggregate all samples from the posterior distributions $p\left({\mu }_j|z_{\mathrm{1,}j}\mathrm{,}\cdots \mathrm{,}{z}_{m\mathrm{,}j}\right)\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}k$, as proposed earlier, rather than using their means.

Figure 2(a) and (b) present the results of the first and second schemes, respectively, with m=200, k=200, N_con=100. We can see that although the frequency distributions obtained by these two schemes have almost the same expectation, their variances are quite different; that is, the second scheme will obtain a variance much closer to the real case. We believe this was caused by the loss of observation information during the averaging step.

Fig. 2Full-screen View

(Color online) Finite cylindrical reactor problem: comparison of schemes 1 and 2, in which scheme 1 aggregates k number of means whereas scheme 2 aggregates k× N number of samples. Not only is the result of scheme 2 intuitively better, but scheme 2 also obtains expectation and variance, particularly the variance, closer to the real values.

More importantly, the estimation of the second scheme is intuitively much better than that of the first. This is because more samples are obtained in the second scheme, resulting in a more robust estimation.

4.1.2

Convergence Analysis

In this subsection, we discuss the effect of the parameters of the method on the result (i.e., the number of sensors m, number of samples k, and number of convergence steps of the Metropolis-Hastings algorithm N_con) and analyze the convergence of the method by setting $N=N_{\text{tot}}-N_{\text{con}}=100$.

Figure 3 shows the effect of N_con on the performance of the proposed method under different settings of (m,k). This shows that the Metropolis-Hastings sampling algorithm used in the ensemble Bayesian method converges rapidly. When N_con is somewhat large, such as N_con=100, the increase hardly affects the performance of the method.

Fig. 3Full-screen View

(Color online) Finite cylindrical reactor problem: the effect of increasing N_con on the results under three sets of parameters (m,k), i.e., (100,100), (150,150), (200,200) respectively. The curves in all three cases show that the effect of N_con on the results is minimal.

Figure 4 shows the effect of k on the performance of the method under different settings of m when N_con = 100 is fixed. This shows that when k is small, there is a fluctuation in performance, and as k increases, the performance tends to stabilize; that is, the expectation and variance tend to be closer to the real values, which are much closer than the prior values.

Fig. 4Full-screen View

(Color online) Finite cylindrical reactor problem: when N_con=100 is fixed, the effect of increasing number of samples k on the results under three sets of parameters m, i.e., 100, 150, 200 respectively. The curves in all three cases show that the performance of the method improves with the increase of m and quickly converges to the real situation.

Figure 5 shows the effect of m on the performance of the method when N_con=100 and k=200 are fixed. This shows that as m increases, the performance of the method improves; that is, the expectation and variance converge to the neighborhoods of the real values, and this expectation trend is particularly pronounced.

Fig. 5Full-screen View

(Color online) Finite cylindrical reactor problem: when $N_{\text{con}}=\mathrm{100,}k=200$ is fixed, the effect of increasing number of sensors m on the results. The curve clearly shows the convergence process of algorithm performance. Compared with N_con and k, m has the greatest influence on the method.

After applying the ensemble Bayesian method by setting $m=\mathrm{200,}k=\mathrm{200,}{N}_{\text{con}}=100$, we obtained the frequency distribution shown in Fig. 2(b). Intuitively, this frequency distribution provides an effective estimation of the real distribution. Equation (19) describes the KL divergence calculations. From the perspective of the KL divergence, we again conclude that the distribution for the estimation derived using the ensemble Bayesian method is much closer to the real distribution than the prior one. $KL\left(p_{\text{prior}}\left|\right|p_{\text{real}}\right)=\mathrm{1.418,}KL\left(p_{\text{est}}\left|\right|p_{\text{real}}\right)=\mathrm{0.180.}$ (19)

4.1.3

Estimation of several parameters

In this section, we test the effectiveness of the method for estimating several parameters simultaneously because several parameters are simultaneously uncertain in many cases. Now, we use the ensemble Bayesian method proposed in this study to simultaneously estimate the distributions of three parameters A,R_e and H_e. Consider the following system: $\{\begin{array}{l}y\left(x\mathrm{,}{\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)={\mu }^{\left(3\right)}\cdot J_0\left(\frac{2.405x}{{\mu }^{\left(1\right)}}\right)\cos \left(\frac{\pi x}{{\mu }^{\left(2\right)}}\right)\hfill \\ z\left(x\mathrm{,}{\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)=y\left(x\mathrm{,}{\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)\left(1+ϵ\right)\hfill \end{array},$ (20) where $x\in [\mathrm{0.5,7.5}]\subset R$ is a one-dimensional spatial position; ${\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}$ are parameters that we are interested in and assumed to have a real joint distribution; J₀ is the zeroth-order first-order Bessel function. $y\in R$ is an observable physical quantity, and $z\in R$ is the observation of y with noise $ϵ\sim N\left(\mathrm{0,}\sqrt{0.2}\right)$.

We experimented with the simulation datasets. To prepare the observational data for our experiments, we first selected 200 positions at equal intervals in [0.5,7.5] and repeatedly obtained observations for the corresponding positions under 200 pairs of (μ⁽¹⁾,μ⁽²⁾,μ⁽³⁾) sampled from the real joint distribution.

Suppose the real joint distribution of (μ⁽¹⁾,μ⁽²⁾,μ⁽³⁾) be 3-dimensional Gaussian distribution with expectation $E\left({\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)$ and covariance matrix $Cov\left({\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)$ as follows: $\begin{array}{l}E\left({\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)=\left(\begin{array}{c}11\\ 18\\ 20\end{array}\right),\\ Cov\left({\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)=\left(\begin{array}{ccc}0.25& 0.1& 0\\ 0.1& 0.25& 0\\ 0& 0& 1\end{array}\right).\end{array}$ It should be noted that the correlations between them are assumed here.

Prior distribution $\pi \left({\mu }_1\mathrm{,}{\mu }_2\mathrm{,}{\mu }_3\right)$ is also a 3-dimensional Gaussian distribution with expectation $E_p\left({\mu }_1\mathrm{,}{\mu }_2\mathrm{,}{\mu }_3\right)$ and covariance matrix $Co{v}_p\left({\mu }_1\mathrm{,}{\mu }_2\mathrm{,}{\mu }_3\right)$ as follows: $\begin{array}{l}{E}_p\left({\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)=\left(\begin{array}{c}12\\ 19\\ 22\end{array}\right),\\ Co{v}_p\left({\mu }^{\left(1\right)}\mathrm{,}{\mu }^{\left(2\right)}\mathrm{,}{\mu }^{\left(3\right)}\right)=\left(\begin{array}{ccc}0.16& 0& 0\\ 0& 0.16& 0\\ 0& 0& 1.44\end{array}\right).\end{array}$

After applying the ensemble Bayesian method by setting m=200, k=200, N_con=100, we obtained three frequency distributions, as shown in Fig. 6. We find that the expectations and variances of all these frequency distributions are significantly close to the expectations and variances of the real marginal distributions. Intuitively, the frequency distributions provide an effective estimation; therefore, the ensemble Bayesian method is also suitable for the simultaneous estimation of several parameters. From the perspective of the KL divergence, the estimation also shows a significant improvement compared to the prior distribution. Equation (21) provides the calculations for the KL divergences. $\begin{array}{l}KL\left(p_{\text{prior}}\left({\mu }_1\right)\left|\right|p_{\text{real}}\left({\mu }_1\right)\right)=2.431\hfill \\ KL\left(p_{\text{est}}\left({\mu }_1\right)\left|\right|p_{\text{real}}\left({\mu }_1\right)\right)=0.300\hfill \\ KL\left(p_{\text{prior}}\left({\mu }_2\right)\left|\right|p_{\text{real}}\left({\mu }_2\right)\right)=2.431\hfill \\ KL\left(p_{\text{est}}\left({\mu }_2\right)\left|\right|p_{\text{real}}\left({\mu }_2\right)\right)=0.130\hfill \\ KL\left(p_{\text{prior}}\left({\mu }_3\right)\left|\right|p_{\text{real}}\left({\mu }_3\right)\right)=2.038\hfill \\ KL\left(p_{\text{est}}\left({\mu }_3\right)\left|\right|p_{\text{real}}\left({\mu }_3\right)\right)=0.039\hfill \end{array}.$ (21)

Fig. 6Full-screen View

(Color online) Finite cylindrical reactor problem: Result of three marginal distributions when setting the number of sensors m=200, number of samples k=200, and number of convergence steps of Metropolis-Hastings algorithm N_con=100.

4.2

the 2D IAEA benchmark problem

[42] gives the definition of the 2D IAEA benchmark problem and [44] gives its implementation with neutronic code. Figure 7 illustrates the geometry of the reactor. Note that only one-quarter is given because the rest can be inferred by symmetry along the x and y axes. We denote this quarter by Ω, which is composed of four subregions with different physical properties: Neumann boundary conditions for the left and bottom boundaries and the mixed boundary condition for the external border.

Fig. 7Full-screen View

(Color online) Geometry of 2D IBP, upper octant: region assignments, lower octant: fuel assembly identification [44].

The 2D IAEA benchmark problem as modeled by two-dimension two-group diffusion equations. Specifically, the flux $\varphi ={\left({\varphi }_1\mathrm{,}{\varphi }_2\right)}^{\text{T}}$ (index 1 and 2 denote the high and thermal energy, respectively) satisfies the following eigenvalue problem [5]: Find $\left(\lambda \mathrm{,}\varphi \right)\in C\times L^{\infty }\left(\text{Ω}\right)\times L^{\infty }\left(\text{Ω}\right)$, s.t. $\{\begin{array}{l}-\nabla \left(D_1\nabla {\varphi }_1\right)+\left({\displaystyle \sum {}_{\text{a}\mathrm{,1}}}+{\displaystyle \sum {}_{1\to 2}}\right){\varphi }_1=\lambda {\chi }_1\left(\nu {\displaystyle \sum {}_{\text{f}\mathrm{,1}}}{\varphi }_1+\nu {\displaystyle \sum {}_{\text{f},2}}{\varphi }_2\right)\hfill \\ -\nabla \left(D_2\nabla {\varphi }_2\right)+{\displaystyle \sum {}_{\text{a}\mathrm{,2}}}{\varphi }_2-{\displaystyle \sum {}_{1\to 2}}{\varphi }_1=\lambda {\chi }_2\left(\nu {\displaystyle \sum {}_{f\mathrm{,1}}}{\varphi }_1+\nu {\displaystyle \sum {}_{\text{f},2}}{\varphi }_2\right)\hfill \end{array},$ (22) where ${\displaystyle \sum {}_{1\to 2}}$ is the macroscopic scattering cross section from groups 1 to 2; Di, ${\displaystyle \sum {}_{\text{a},i}}$, ${\displaystyle \sum {}_{\text{f},i}}$, and χi are the diffusion coefficient, macroscopic absorption cross section, macroscopic fission cross section, and fission spectrum of group i, $i\in \left\{\mathrm{1,2}\right\}$, respectively; ν is the average number of neutrons emitted per fission. In other words, these quantities are the model parameters. The general settings are listed in Table 1.

Table 1Full-screen View

General settings of the parameters of the IAEA 2D benchmark problem.

Region	D1 (cm)	D2 (cm)	${\displaystyle \sum {}_{1\to 2}}$ (cm-1)	${\displaystyle \sum {}_{\text{a},1}}$ (cm-1)	${\displaystyle \sum {}_{\text{a},2}}$ (cm-1)	ν ∑f,1 (cm-1)	ν∑f,2 (cm-1)	χ1	χ2	Material
Ω1	1.50	0.40	0.02	0.01	0.080	0.00	0.135	1	0	Fuel 1
Ω2	1.50	0.40	0.02	0.01	0.085	0.00	0.135	1	0	Fuel 2
Ω3	1.50	0.40	0.02	0.01	0.130	0.00	0.135	1	0	Fuel 2 + Rod
Ω4	2.00	0.30	0.04	0.00	0.010	0.00	0.000	0	0	Reflector

Show more

Later, we assumed D₂ in Ω₄ was uncertain and used the proposed ensemble Bayesian method to estimate its distribution, as it is probably one of the most important parameters. We assume that all other parameters are certain and set them according to Table 1 except for D₂ in Ω₄. We denote this as μ, which is what we are interested in.

We experimented with simulation datasets. To prepare the observation data for our experiments, we first selected 45 positions (i.e., (10,10),(30,10),(30,30),…, (170,10),…,(170,170)) in area Ω and repeatedly obtained the observations for the corresponding positions under 200 pairs of μ sampled from the real distribution.

Moreover, the two-dimensional two-group diffusion Eq. (22) must be solved in the algorithm process. This was achieved by employing the generic high-quality finite element solver FreeFem++ [45].

Suppose that the real distribution of μ is a Gaussian distribution with expectation $E(\mu )$ and variance $Var(\mu )$ as follows: $E\left(\mu \right)=\mathrm{2.00,}Var\left(\mu \right)=\mathrm{0.04,}$ (23) where $E(\mu )$ is the value listed in Table 1.

The prior distribution π (μ) is also a Gaussian distribution with expectation $E_{\text{p}}\left(\mu \right)$ and variance $Va{r}_{\text{p}}\left(\mu \right)$ as follows: $E_{\text{p}}\left(\mu \right)=\mathrm{2.20,}Va{r}_{\text{p}}\left(\mu \right)=\mathrm{0.05.}$ (24)

After applying the ensemble Bayesian method with setting N_con=100, we obtained the frequency distribution shown in Fig. 8. We find that the expectation of the frequency distribution is significantly close to the real expectation, and the variance is also closer to the real variance than the prior distribution. Furthermore, the frequency distribution intuitively provided an effective estimation. From the perspective of the KL divergence, the estimation also shows a significant improvement compared to the prior estimation. Equation (25) describes the KL divergence calculations. In this case, we used the engineering software Freefem++ to implement our algorithm because this model only has numerical solutions; therefore, the proposed ensemble Bayesian method has potential for engineering applications. $KL\left(p_{\text{prior}}\left|\right|p_{\text{real}}\right)=\mathrm{0.5134,}KL\left(p_{\text{est}}\left|\right|p_{\text{real}}\right)=\mathrm{0.0213.}$ (25) For the finite cylindrical reactor, the computation time of one dimensional case (Fig. 2(b)) is 111.7 s; the computation time of three dimensional case (Fig. 6) is 359.6s. For the 2D IAEA benchmark (Fig. 8), the computation time is about 130000s. The main reason why 2D IAEA benchmark takes more time is that FreeFem++ takes up most of the runtime, and the most important problem is that FreeFem++ is not called or connected to the code efficiently enough. Therefore, if there is a way to call FreeFem++ efficiently, it can significantly save computing time.

Fig. 8Full-screen View

(Color online) 2D IAEA benchmark problem: Result of the method when setting the number of samples k=200 and the number of convergence steps of Metropolis-Hastings algorithm N_con=100.

Conclusion

In this study, we considered model parameters from the perspective of random variables in the context of nuclear reactor engineering and proposed a general form of parameter distribution inference problem. In the context of this parameter distribution estimation problem, we conducted a preliminary exploration and proposed an ensemble Bayesian method to estimate the parameters by obtaining frequency distributions, combined with a special type of experimental setup. Simultaneously, we introduced the KL divergence theory to quantify the estimation performance of the method. Various numerical experiments were conducted, including a finite cylindrical reactor model, which has an analytical solution, and the 2D IAEA benchmark problem, which only has a numerical solution. From the results of these tests, the following conclusions were drawn:

• For those two models, the frequency distributions intuitively provide effective estimation for the parameters. The expectations are significantly close to the real ones, and the variances are also more accurate than the prior distributions. From the perspective of KL divergence, the method also has good performance;

• The ensemble Bayesian method can estimate several parameters simultaenously, even if there are correlations between them;

• When the model has only numerical solution, we use engineering software Freefem++ to implement our algorithm. The method also works well in this case, which means it has potential for engineering applications;

• The convergence speed of the method is fast. Thus, in general, the ensemble Bayesian method we propose has optimistic application prospects.

The proposed method shows high potential for engineering applications to correct prior distributions when using the data assimilation technique for parameter estimation.

Further studies could investigate the performance of the proposed method under moderate-scale parameters, the performance for non-Gaussian parameters, and characteristic parameters that are likely to be difficult to reproduce. As an important review of machine learning (ML) in nuclear physics, [46] not only describes the different methodologies used in ML algorithms and techniques and some applications, particularly for low- and intermediate-energy nuclear physics, but also provides a valuable summary and outlook on the possible application directions and improvements of ML algorithms in low- and intermediate-energy nuclear physics. Inspired by the review, we will improve the efficiency of the method, and our initial idea is to improve the efficiency of the Metropolis-Hastings sampling method used in this study. Furthermore, we will determine how to involve more physics-informed or physics-guided to develop new machine learning algorithms.