Unfolding neutron spectra from water-pumping-injection multilayered concentric sphere neutron spectrometer using self-adaptive differential evolution algorithm

NUCLEAR ELECTRONICS AND INSTRUMENTATION

Unfolding neutron spectra from water-pumping-injection multilayered concentric sphere neutron spectrometer using self-adaptive differential evolution algorithm

Rui Li，

Jian-Bo Yang ，

Xian-Guo Tuo，

Jie Xu，

Rui Shi

Nuclear Science and Techniques

Vol.32, No.3

Article number 26

Published in print 01 Mar 2021

Available online 16 Mar 2021

DOI：10.1007/s41365-021-00864-4

47701

A self-adaptive differential evolution neutron spectrum unfolding algorithm (SDENUA) is established in this study to unfold the neutron spectra obtained from a water-pumping-injection multilayered concentric sphere neutron spectrometer (WMNS). Specifically, the neutron fluence bounds are estimated to accelerate the algorithm convergence, and the minimum error between the optimal solution and input neutron counts with relative uncertainties is limited to 10^-6 to avoid unnecessary calculations. Furthermore, the crossover probability and scaling factor are self-adaptively controlled. FLUKA Monte Carlo is used to simulate the readings of the WMNS under (1) a spectrum of Cf-252 and (2) its spectrum after being moderated, (3) a spectrum used for boron neutron capture therapy, and (4) a reactor spectrum. Subsequently, the measured neutron counts is unfolded using the SDENUA. The uncertainties of the measured neutron count and the response matrix are considered in the SDENUA, which does not require complex parameter tuning or an a priori default spectrum. The results indicate that the solutions of the SDENUA agree better with the IAEA spectra than those of MAXED and GRAVEL in UMG 3.1, and the errors of the final results calculated using the SDENUA are less than 12%. The established SDENUA can be used to unfold spectra from the WMNS.

Water-pumping-injection multilayered spectrometerNeutron spectrum unfoldingDifferential evolution algorithmSelf-Adaptive control

1 Introduction

Since the first introduction of the Bonner sphere spectrometer (BSS) in 1960 [1], it has been widely used in neutron spectrometry measurements, e.g., the isotopic neutron source [2], boron neutron capture therapy (BNCT) [3-4], and radiation protection near a reactor [5-7], owing to its advantages in isotropic response and wide energy range. A newly designed neutron spectrometer, i.e., the water-pumping-injection multilayered concentric sphere neutron spectrometer (WMNS), uses water as a moderator [8-10], and its principle of neutron spectrometry measurement is similar to that of the BSS. The structure of the WMNS is illustrated in Fig. 1. Seven stainless steel spherical shells are arranged concentrically to build six spherical gaps, five of which are used to contain water (the thickness of the water gap from the outside to the inside is 2.5, 3.75, 8, 1.25, and 1 cm, in that order), and the remaining gap is filled with lead. The water is independently pumped into each gap to form a measurement unit (combination) to moderate the incident neutrons, and an easy-to-replace spherical ³He proportional counter (model: LND 2705) is placed in the innermost part to detect thermal neutrons [11-12]. Five gaps afford up to 32 measurement combinations with different water thicknesses (similar to the ball with different diameters in the BSS); therefore, 32 measurement combinations can be used to obtain the neutron count. Switching between measurement combinations was realized using an external water-pumping-injection system [13]. A 1 cm lead was utilized to measure high-energy neutrons. The WMNS is a portable and flexible neutron spectrometer. An active or passive detector can be used depending on the measurement environment, and only one reading electric system is required. Moreover, the “ball–ball interference” problem in the traditional BSS is eliminated by integrating the multiconcentric spheres. The measurement combinations are coded by 0 and 1, e.g., the combination code of Fig. 1 is “00Pb110,” where “0” represents a gap filled with air, “1” represents a gap filled with water, and “Pb” represents lead.

Fig. 1

Schematic diagram of WMNS

The readings of the ³He proportional counter, also known as measured counts, are the nuclear reaction event counts of ³He(n, p)³H under different measurement combinations. The target spectrum is unfolded from these measured counts using the neutron unfolding algorithm, and the neutron unfolding process can be presented in a discrete form as follows [14]:

C_{j}^{meas} + ε_{j} = \sum_{i = 1}^{n} R_{i j} φ_{i} j = 1, 2, 3, ... m,

(1)

where $C_{j}^{meas}$ is the measured neutron count reading from the jth measurement combination, $ε_{j}$ the reading uncertainty of the jth measurement combination, $R_{i j}$ the response of the jth measurement combination to the neutron of the ith energy group, and $φ_{i}$ the neutron fluence of the ith energy group. Typically, the number of measurement combinations is significantly smaller than the number of energy groups.

In the WMNS, to minimize the time required for switching measurement combinations, 18 measurement combinations were selected. Thirty-six energy groups were divided logarithmically at equal intervals in the range from 10⁻⁹ MeV to 20 MeV to reduce the underdetermined degree of the unfolding process. Therefore, m = 18 and n = 36. The FLUKA Monte Carlo code [15] was used to calculate the response matrix, as shown in Fig. 2, where parallel monoenergetic neutron beams, starting from a disk with a diameter of 28 cm, which was the same as the diameter of the outermost stainless-steel shell of the WMNS, irradiated the WMNS to obtain responses. The distance between the source and the center of the spectrometer was 60 cm. The RESNUCLEi card was used to score the stopping nuclei on ³He. The stopping nuclei were tritium nuclei and protons because of the ³He(n, p)³H reaction; subsequently, half of all the nuclei were collected as the reading of the detector because when each tritium nucleus or proton was produced, a neutron was detected simultaneously. The thermal neutron scattering data S(α, β) were applied to the transport of neutrons below 4 eV in polyethylene and water using the LOW-NEUT and LOW-MAT cards [16].

Fig. 2

Response matrix of WMNS

Currently, various methods, such as the maximum entropy method [17] and iterative method [18] used in the unfolding with Maxed and Gravel 3.1 (UMG 3.1) [19] can be used to unfold the neutron spectrum. Their solutions are compared with those of the present study herein. An excellent priori default spectrum is required when using UMG 3.1 to unfold the neutron energy spectra [17] because in the maximum entropy method, a priori default spectrum is a benchmark for UMG 3.1 and determines when to output the solution; in the iterative method, the priori default spectrum is the initial of the iteration. Shahabinejad et al. [20] used a two-step genetic algorithm (TGA) to unfold neutron energy spectra, and the results showed a closer match in all energy regions, particularly in the high-energy regions, compared with the typical genetic algorithm (GA). Energy groups in the high energy range were unfolded at the first step and then used to construct the initial value of the second step. One year later, they used a particle swarm optimization algorithm (PSOA) [21] to unfold the neutron spectrum from a pulse height distribution and a response matrix, and the results demonstrated consistency with the TGA. In the PSOA, the acceleration constants c₁, c₂, and inertia weight w are empirically predefined by the authors. Hoang et al. [22] applied a different two-step GA to unfold neutron spectra obtained from activation foils. Unlike a previous study [20], in the first step, only the energy groups in the region from 20 to 35 MeV were unfolded, whereas in the second step, the entire energy spectrum was unfolded while maintaining the result of the first step. K. Chang et al. [9] established a backpropagation artificial neural network [23] neutron spectrum unfolding code, and the neural network was trained under 32 neutron spectra. Furthermore, the capability of the code was verified using eight neutron spectra. As mentioned above, the methods in UMG 3.1 rely on an a priori default spectrum, whereas parameter tuning in GA and PSOA frameworks are complicated, and neural network training is time consuming and complex.

In this study, we focused on the necessity to unfold the neutron spectrum from a WMNS. A self-adaptive differential evolution neutron spectrum unfolding algorithm (SDENUA), which includes the neutron fluence bound estimation and parameter self-adaptive control technique, is proposed herein. The error between the input neutron counts and the calculated counts was limited to 10^-6 to improve the quality of the solutions and reduce the calculation time. The measured neutron counts of 1) the spectrum of Cf-252 and 2) its spectrum after being moderated, 3) a spectrum used for BNCT, and 4) a spectrum from a reactor in the IAEA 403 report [24] were simulated using the FLUKA code. Subsequently, the measured neutron counts with relative uncertainties were obtained using the SDENUA. The SDENUA does not require complex parameter tuning and an a priori default spectrum, and the established SDENUA can be used unfold spectra from the WMNS.

The remainder of this paper is organized as follows. In the “Material and methods” section, first, the neutron fluence bounds are estimated; second, the techniques related to the self-adaptive differential evolution algorithm used for neutron energy spectra unfolding are explained in detail along with the formulation of the algorithm; and third, the termination criterion of the algorithm is proposed. In the “Results and discussion” section, the unfolded spectra of this study are shown and compared with those of the UMG 3.1, and the uncertainties of the unfolded spectra are discussed. Finally, the conclusion is presented in the “Conclusion” section.

2 Material and methods

2.1 Estimation of neutron fluence bounds

In the optimization problem, the search space bound of each variable must be pre-estimated because the scale of the search space significantly affects the operating time and convergence of the algorithm [14]. In other words, the neutron fluence of each energy group must be bounded before unfolding. It is well known that the actual neutron energy spectrum is due to its physical properties, and that the minimum neutron fluence value of all energy groups is 0; however, the upper fluence must be estimated. In Ref. [25], it was assumed that the measured count of a particular measurement unit is fully contributed by a particular energy group, and that the minimum fluence value of these estimates is the upper bound of the energy group, as shown in Eq. (3). For the ideal monochromatic pulse neutron energy spectrum, this method can directly provide the neutron fluence at which the pulse is located. Although this method is based on strict mathematical derivation, because the contribution of the fluence outside the particular energy group to the neutron count is disregarded, the result of this method is only approximate.

The neutron energy spectrum is typically continuous [14], and we assume that the fluence change between adjacent energy groups is relatively smooth. Therefore, the range estimated using the method above is narrowed by the following conditions:

φ_{i}^{\max} = \frac{φ_{i}^{1} + φ_{i}^{3}}{2},

(2)

φ_{i}^{1} = M I N (C_{j}^{meas} / R_{i j}),

(3)

φ_{i}^{3} = \frac{φ_{i}^{1} \cdot φ_{⌈ i / 3 ⌉}^{1}}{\sum_{i}^{i + 2} φ_{i}^{1}},

(4)

where $φ_{i}^{\max}$ is the fluence upper bound of the ith energy group, $φ_{i}^{1}$ the fluence upper bound of the ith energy group estimated using the method in Ref. [25], and $φ_{i}^{3}$ the fluence upper bound of the ith energy group with three energy groups; $⌈ i / 3 ⌉$ rounds to the nearest integer greater than or equal to i/3, which is the first of every three $φ^{1} = [φ_{1}^{1}, φ_{2}^{1}, ... φ_{i}^{1} ... φ_{n}^{1}]$ (n is the number of the energy group), whereas $C_{j}^{meas}$ and $R_{i j}$ are the same as those in Eq. (1). Finally, the fluence upper bound of the ith energy group is in the interval $φ_{i} \in (0, φ_{i}^{\max})$ .

It is noteworthy that a smooth neutron spectrum is a prerequisite; otherwise, the fluence of the peak will be underestimated. Therefore, estimating the fluence upper bound of the energy spectrum, which contains a sharp peak, is not recommended.

2.2 Self-adaptive differential evolution algorithm

The differential evolution algorithm was introduced by Storn et al. in 1996 [26]. This algorithm has garnered significant attention owing to its simple framework and robust global search capabilities, and the agreement between the individuals and solutions is evaluated by the fitness value generated by the fitness function. The iterative loop of the algorithm includes evolutionary operations, such as initialization, mutation, crossover, and selection. For neutron unfolding, the algorithm population is composed of several individuals. An individual corresponds to a neutron spectrum comprising several genes. The positions of the genes correspond to the positions of the energy group, and the values of the genes correspond to the neutron fluence, which is the variable of the unfolding problem to be solved. The SDENUA is introduced in detail below.

Initialization of the neutron fluence of each energy group was achieved by randomly selecting the neutron fluence from the neutron fluence estimation interval $(0, φ_{i}^{\max})$ . A key issue in initialization is the number of individuals included in the population. An excessively small population size can result in premature convergence, whereas an excessively large population size will result in a long calculation time and insufficient mixing. Gämperle et al. analyzed the results of a 20-dimensional (20 variables) problem [32] and reported that a reasonable option for the population size was between three and eight times the number of variables; meanwhile, 10 times the number of variables was recommended in Ref. [27]. To ensure that the population contained a sufficient number of diverse vectors to participate in the mutation, as well as to improve the search and traversal ability of the population in the evolutionary process, the population size in this study was set to 10 times the number of the energy group (PopSize = 10 × 36), and the population size was maintained throughout the entire unfolding process. Such a population size implies a larger number of calculations; hence, a time-conserving technique is proposed, as will be described in Sect. 2.3.

The mutation operation, which guides the direction of evolution of the population, was executed using Eq. (5) [28]. The search step length is controlled by the scaling factor. Therefore, the mutation operation provides two functions: the search direction and search step length control evolution.

v_{i, g} = k_{i, g} + F_{i} \cdot (x_{best, g}^{B * P o p S i z e} - k_{i, g}) + F_{i} \cdot (x_{r, g} - x_{f, g}),

(5)

where $v_{i, g}$ is the ith temporary individual in the gth generation, $k_{i, g}$ the ith target individual in the gth generation, and $x_{best, g}^{B * P o p S i z e}$ the high-fitness individual. $x_{r, g}$ is randomly selected from the current population P, $x_{f, g}$ is randomly selected from ( $P_{f} \cup P$ ), and $P_{f}$ is a set used to save failed individuals (i.e., individuals with lower fitness in the selection step). $F_{i}$ is the scaling factor for each target individual.

Prior to the mutation operation, all individuals in the current population P were sorted in the descending order based on their fitness. High-fitness individuals were randomly selected from the top 100*B% individuals after sorting, and B was a uniform random number on the interval [0.05, 0.60]. The upper bound at 0.6 of B was more suitable for the 36-dimension unfolding problem as it enabled more high-fitness individuals to participate in the mutation step as well as reduced the risk of premature convergence. The set $P_{f}$ was created with a size of 100 from the beginning of the first iteration, and the individuals with lower fitness in the selection step were sorted as follows: When $P_{f}$ was full, the failed individuals in $P_{f}$ were updated based on the “first in, first out” rule. $P_{f}$ maintained the failed individuals from 100 generations and provided more different genes for mutation.

Each scaling factor was generated independently based on a normal distribution, i.e., $F_{i} ~ N (u_{F}, 0.1)$ , which can yield better results compared with the Cauchy sampling method, and $u_{F}$ was updated [28] at the end of each generation as follows:

u_{F} = (1 - c) \cdot u_{F} + c \cdot {mean}_{L} (S_{F}),

(6)

where mean_L(·) is the Lehmer mean, expressed as

{mean}_{L} (S_{F}) = \frac{\sum_{F \in S_{F}} F^{2}}{\sum_{F \in S_{F}} F} ，

(7)

where $S_{F}$ is a set for all failed mutation factors in each generation, and the $u_{F}$ saved in $S_{F}$ is used to provide the prior historical information for the next $u_{F}$ update. $S_{F}$ is blanked at the beginning of each generation. c is a uniform random number on the interval [0.05, 0.20], and it controls the life span of $u_{F}$ with generations from 5 to 20 [28]; $u_{F} =0 .5$ at the initialization.

The crossover operation is based on the temporary individuals generated by the mutation operation and causes them to crossover with the target individuals as follows:

u_{i, g}^{j} = {\begin{cases} v_{i, g}^{j}, r^{j} \leq C R_{i} \\ k_{i, g}^{j}, otherwise \end{cases},

(8)

where $u_{i, g}^{j}$ is the jth gene of the ith candidate individual, $v_{i, g}^{j}$ the jth gene of $v_{i, g}$ , $k_{i, g}^{j}$ the jth gene of $k_{i, g}$ , $r^{j}$ a uniform random number on the interval [0, 1], and $C R_{i}$ the crossover probability of each candidate individual.

The basic unit of crossover is the gene, and the genes in the temporary and target individuals are extracted to construct a candidate individual based on the crossover probability. The candidate individuals indirectly transmit the effect from the mutation operation to the target individual, and the effect is controlled by the crossover probability. The self-adaptive adjustment of the crossover probability based on historical evolution information is more effective than the traditional constant control method.

Each $C R_{i}$ is generated, as shown in Eq. (9), and updated using Eq. (10) [28].

C R_{i} = {\begin{cases} C R_{i} ~ N (u_{C R}, 0.1), otherwise \\ C R_{i} = 0.5 ， C R_{i} < 0.5 \\ C R_{i} = 0.95 ， C R_{i} > 0.95 \end{cases},

(9)

u_{C R} = (1 - c) \cdot u_{C R} + c \cdot S_{C R}^{M e a n},

(10)

where $S_{C R}^{Mean}$ is the average of all elements in the set, $S_{C R}$ , and each $C R_{i}$ is sent to the set $S_{C R}$ if the fitness of the target individual is higher than that of the candidate individual in the selection operation. $S_{C R}$ is blanked at the beginning of each iteration. c is a uniform random number on the interval [0.05, 0.20] and controls the life span of $u_{F}$ with generations from 5 to 20 [28]; at the initialization [26, 28]. To ensure that the genes of both the temporary and target individuals can be passed on partially to the candidate individuals, $C R_{i}$ was truncated to [0.5, 0.95] [29, 30]. The candidate individual generated by the crossover operation was then sent to the selection operation.

The selection operation determines whether to refuse or allow the candidate individual to enter the population based on the following:

k_{i + 1, g} = {\begin{cases} u_{i, g}, f (u_{i, g}) > f (k_{i, g}) \\ k_{i, g}, otherwise \end{cases},

(11)

where f (·) is the fitness function. If the candidate individuals with lower fitness are rejected from the population to maintain the average fitness of the population at a higher level. In that case, the mutation and crossover operations will be futile if we do not utilize the candidate individuals.

The fitness function [22] is used to evaluate the closeness between the input neutron counts and individuals (solutions), as follows:

f = {[\sum_{j = 1}^{m} \frac{{(C_{j}^{u} - C_{j}^{cal})}^{2}}{{(C_{j}^{u})}^{2}}]}^{-1},

(12)

where $C_{j}^{u}$ is the input neutron count of the jth measurement combination generated with a relative uncertainty, $C_{j}^{cal} = \sum_{i}^{n} R_{i j}^{u} φ_{i}^{cal}$ is the calculated neutron count of the jth measurement combination, $R_{i j}^{u}$ is a response generated with relative uncertainty, and $φ_{i}^{cal}$ is the neutron fluence of the ith energy group of the calculated spectrum (solution). If the solution is closer to $C_{j}^{u}$ , then the fitness is higher. In the spectrum unfolding process, the fitness function is the only criterion for determining the quality of the solution. Therefore, the performance of the fitness functions can significantly affect the solution as well as the performance of the algorithm. To reduce the dependence of the algorithm on the fitness function, the fitness value of the optimal solution was restricted to obtain a physically acceptable solution, as will be described in the next section.

2.3 Termination criterion

The neutron energy spectrum unfolding problem is a first-kind Fredholm integral problem, and a perfect solution cannot be obtained based on the integral fitness function shown in Eq. (12); furthermore, a fundamental hypothesis for the unfolding algorithm is that solutions with an acceptable spectral quality can be obtained based on the fitness function [22]. In addition, to improve the probability of convergence and obtain high-fitness solutions, a larger maximum iteration number is typically required to truncate the fitness of the final solution, such as that presented in [14, 20], where the authors empirically defined a larger maximum iteration number for spectrum unfolding. However, an overestimation of the maximum iteration number would increase the calculation time, whereas an underestimation would yield a pseudo-optimal solution before convergence. Therefore, a better termination criterion is proposed herein.

The spectrum quality factor (QS) [14] was used as a metric to evaluate the quality of the solution; it is expressed as

Q S =100 \cdot \sqrt{\frac{\sum_{i = 1}^{n} {(φ_{i}^{ref} - φ_{i}^{cal})}^{2}}{\sum_{i = 1}^{n} {(φ_{i}^{cal})}^{2}}},

(13)

where $φ_{i}^{ref}$ is the neutron fluence of the ith energy group of the reference spectrum (solution) and $φ_{i}^{cal}$ is the same as that in Eq. (12). A perfect solution results in a QS value of 0.

In the unfolding process, the quality of the optimal solution is expected to improve as the fitness of the solutions increases; this implies that the QS of the solutions conforms to the monotonic non-increasing trend. The relationship between QS and the fitness of the “for BNCT” spectrum is shown in Fig. 3; as shown, when the fitness of the solutions increased, the QS of the solutions fluctuated. Although the fitness of the final solution was extremely high, the QS of the final solution was not the lowest in history. In other words, the optimal solution was disregarded, and the calculation resulted in a QS rebound after the generation with the lowest QS, which was meaningless and harmful in some cases. However, the QS is an evaluation indicator based on a known energy spectrum. In the actual energy spectrum unfolding process, we can only use fitness to evaluate the solutions. To minimize the QS of the final solution, after numerous experiments, the fitness of the optimal solution was limited to 10⁶ as the iteration termination criterion to replace the maximum iteration number; this reduced the operating time and improved the quality of the final optimal solution.

Fig. 3

Relationship between QS and fitness of “for BNCT” spectrum

To investigate the effects of the uncertainties from the measured neutron counts and response matrix on the final results, the inputs of each run in the SDENUA were generated as follows:

C_{j}^{u} = C_{j}^{meas} + r a n d (- ε_{j}, ε_{j}) C_{j}^{meas},

(14)

R_{i j}^{u} = R_{i j} + r a n d (- σ_{i j}, σ_{i j}) R_{i j},

(15)

where $C_{j}^{u}$ and $R_{i j}^{u}$ are the measured neutron count with uncertainty estimation and the response matrix with uncertainty estimation, respectively; $C_{j}^{meas}$ , $R_{i j}$ , and $ε_{j}$ are the same as those in Eq. (1), $σ_{i j}$ the uncertainty of each response, and rand(·) the uniform random sampling function. The maximum errors of the measured neutron count simulated from the relative reference spectra were 0.37%, 0.75%, 2.23%, and 3.10%. The average and maximum errors of the response matrix were 0.75% and 10.23%, respectively. The flowchart of the SDENUA is shown in Fig. 4.

Fig. 4

Flowchart of SDENUA

3 Results and discussion

As shown in Fig. 5, after multiple independent runs, the QS of the solutions with an average of 20 independent runs was lower except in the “reactor” spectrum. Therefore, 20 times is a reasonable setting for actual applications.

Fig. 5

QS of solutions with average over multiple independent runs in SDENUA

As shown in Fig. 6, the termination generations of the final optimal solution of each spectrum differed for 20 runs when the fitness of the solutions reached 10⁶. The average termination generation for 20 runs of the “Cf-252” spectrum was the 1157th generation; “Cf-252 Mod” spectrum, the 544^th; “for BNCT” spectrum, the 451^th; and “reactor” spectrum, the 472^th. The termination generations suitable for “Cf-252” and “Cf-252 Mod” were overestimated, whereas those for “for BNCT” and “reactor” were underestimation. Therefore, a universal maximum number of iterations might not exist. Hence, it is proved once again that in the neutron spectrum unfolding problem, particularly when considering multiple energy spectrum types, using the maximum number of iterations as the termination condition is not an optimal solution.

Fig. 6

Termination generation of final optimal solution of 20 runs when fitness of solutions reached 10⁶

The termination generations varied significant between runs, and more generations were required for the Cf-252 spectrum. As shown in Fig. 8a, many energy groups with zero fluences appeared. To ensure a small error between the input and calculated counts, the algorithm can only render the fluences approximately 0 in the side positive range because the negative fluences were rejected. This implies that the solution space that can satisfy the error constraint is narrowed, and more iterations are required. Moreover, only the energy groups with fluences can contribute to the neutron counts when the solution is convolved with the response matrix. The fewer nonzero terms in the Cf-252 solution, the greater is the instability of the Cf-252 solution; hence, the distribution of the solutions will expand.

Fig. 8

Neutron spectra unfolding results (with average of over 20 independent runs in SDENUA): (a) spectrum of isotope source of Cf-252, (b) spectrum of Cf-252 source after being moderated, (c) spectrum used for BNCT, and (d) spectrum of a certain reactor from Germany. Uncertainties for the calculated spectra were less than 5%. (Calculations are presented in the Appendix).

The oscillations in the unfolded spectra are an inherent feature of the numerical solution to the Fredholm equation [19, 20, 30], and a Gaussian smoothing method was adopted to smooth the final optimal solution to overcome the oscillations. In Gaussian smoothing, the smoothing window size is 7, and sigma is 1.4.

As shown in Fig. 8, the neutron fluence at the peak position of the solutions of the SDENUA were underestimated; this is attributable to the optimal solution of the SDENUA smoothed by the Gaussian method increasing the error artificially. However, the results shown in Fig. 7 indicate that the positive effect of Gaussian smoothing was greater than that the smoothing effect.

Fig. 7

QS of final solutions from SDENUA smoothed using Gaussian and without smoothing

As described above in the introduction section, the methods in UMG 3.1 begin with an a priori default spectrum. For a fair comparison, an excellent priori spectrum was obtained using the UMG 3.1 unfolding code. Furthermore, the inputs ( $C_{j}^{u}$ and $R_{i j}^{u}$ ) were generated using Eqs. (14) and (15), respectively, and the chi-square was set to 10⁻⁶ to compare with the fitness upper bound in the SDENUA. Fig. 8 shows the priori spectra and the unfolding results of the four spectra in the IAEA 403 report [24].

Except for Fig. 8a, the unfolding results of the other three energy spectra show obvious errors compared with the reference spectra in the energy range from 10^-9 to 10^-7MeV. As shown in Fig. 2, the response functions of energies from 10^-9 to 10^-7 MeV overlapped significantly, thereby weakening the unfolding power of the response functions in this energy range. Owing to the integral fitness function, as shown in Figs. 8c and 8d, although the final optimal solutions fluctuated in a large range around the reference spectra, the fitness of the final solutions still reached 10⁶. This implies that the solution space that can satisfy the error constraint is expanded, and that fewer iterations for searching are required, as shown in Fig. 6.

The results yielded by UMG 3.1 are interesting. Although an excellent priori default energy spectrum was obtained, the agreement between the results and reference spectra was poor. This may be due to the neutron counts and response matrix input with uncertainties, as well as UMG 3.1 unfolding the neutron spectra based strictly on the inputs.

The uncertainties of the unfolded results were from both the uncertainty terms of the measured neutron counts and the response matrix. Fig. 9 shows that the final results yielded by the SDENUA had errors ranging from -6% to 12%, as calculated using the following equation:

Fig. 9

Uncertainties of unfolded results. Uncertainties for calculated neutron counts were less than 5%. (Calculations are presented in the Appendix).

E_{j} = \frac{C_{j}^{cal} - C_{j}^{meas}}{C_{j}^{meas}} \times 100,

(16)

where $E_{j}$ is the error of the jth measurement combination, $C_{j}^{cal}$ is the calculated count of the jth measurement combination, and $C_{j}^{meas}$ is the same as that in Eq. (1).

Although the fitness value of the final solution was limited to 10⁶, a 12% error existed between the unfolding result and the exact measured count owing to the introduction of uncertainties from both the measured counts and the response matrix. In the neutron unfolding problem, more attention is typically focused on obtaining a suitable spectrum instead of a spectrum with ultra-small counting errors [20]. This is because ultra-small errors do not yield an ideal analytical solution for the first-kind Fredholm integral problem biased on the integral fitness function [14]. Therefore, combining Figs. 8 and 9, we believe that the errors of the unfolded results show an acceptable accuracy level.

The fitness function in this study was used to match the integral quantities of the input neutron counts and the calculated neutron counts. This will trigger a discussion regarding the mechanism of the selection operation. The candidate individual that has failed to evolve will be denied entry into the population in time to avoid a deterioration in the population quality. Meanwhile, as shown in Fig. 3, when the fitness value of an individual (solution) is higher, the corresponding QS (the error indicator) may be low; hence, the error between the solution and the actual energy spectrum may be worse even though the fitness value is high. In addition, an individual is composed of multiple genes, and a few excellent genes will be immediately eliminated from the population because of the low fitness of the candidate individuals. Consequently, a potential candidate individual will be immediately affected because of its low fitness value.

Regarding the SDENUA established in this study, the QS values of the optimal solutions of “Cf-252” and “Cf-252 Mod” in Fig. 7 were worse than those of “for BNCT” and “reactor” when the final optimal solution reached 10⁶. In other words, limiting the fitness of the final optimal solution might jeopardize the unfolding accuracy of the “Cf-252” and “Cf-252 Mod” energy spectra. We hypothesize that one of the reasons for this phenomenon is that the “Cf-252” and “Cf-252 Mod” energy spectra contained fewer energy groups with fluences, i.e., equivalent to reducing the number of effective constraint items for neutron counts calculated in the convolution process. In summary, for such an energy spectrum containing a large number of 0 fluence energy groups, 10⁶ as the fitness upper of the optimal solution may not be sufficiently large; however, considering other types of energy spectra, a tradeoff was necessary.

Regarding the a priori spectrum in the neutron energy spectrum unfolding problem, researchers have posed different opinions. The developers of UMG 3.1 reported [17] that the problem of neutron energy spectrum unfolding should be based on an excellent priori default energy spectrum. They believed that adding the physical information of the neutron energy spectrum enabled more accurate unfolding results to be obtained. However, in some studies, researchers who used artificial intelligence algorithms reported that it was difficult to estimate an excellent priori spectrum in some cases. Therefore, dependence on the a priori energy spectrum should be reduced. The former was explained from the perspective of physics, whereas the latter from mathematics. We believe that it is difficult to obtain an accurate solution by only solving the first-kind Fredholm integral problem, although this has been attempted previously [19-21]. In terms of the algorithm based on the fitness function, some physical information regarding the neutron energy spectrum can be added to the fitness function as a constraint instead of attempting to construct an excellent priori spectrum, such as the continuity of the neutron spectrum. This might contribute positively to unfolding, which require numerous experiments and in-depth investigations.

4 Conclusion

The SDENUA yielded promising results for a WMNS. In the mutation operation, the information of high-fitness and failed individuals was used to improve the guidance for the evolution direction of the population; hence, the fitness of the final optimal solution reached 10⁶ rapidly to output an acceptable solution, thereby improving the quality of the solution and reducing the operating time. Historical experience information was adopted to perform self-adaptive control of the scaling factor and crossover probability. The constructed self-adaptive difference algorithm was used to unfold the readings simulated from the (1) spectrum of Cf-252 and (2) its spectrum after being moderated, (3) a spectrum used for BNCT, and (4) a spectrum from a reactor in the IAEA 403 report. The unfolding spectra agreed better with the reference spectra than those of UMG 3.1. This demonstrated that in the absence of an a priori default spectrum, and with the uncertainties of the measured neutron counts and response matrix, the unfolded results were at an acceptable level with errors less than 12%.

Reference

R.L. Bramblett, R.I. Ewing, T.W. Bonner,

A new type of neutron spectrometer

. Nucl. Instrum. Methods 9, 1-12 (1960). doi: 10.1016/0029-554X(60)90043-4