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

Rui Li· Jian-Bo Yang,,3· Xian-Guo Tuo,,3 · Jie Xu· Rui Shi

Abstract 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 are 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.

Keywordst Water-pumping-injection multilayered spectrometer · Neutron spectrum unfolding · Differential evolution algorithm · Self-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 spherical3He 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.

The readings of the3He proportional counter, also known as measured counts, are the nuclear reaction event counts of3He(n, p)3H 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]:

Fig. 1 Schematic diagram of WMNS

Fig. 2 (Color online) Response matrix of WMNS

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-9to 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 on3He. The stopping nuclei were tritium nuclei and protons because of the3He(n, p)3H 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].

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 highenergy 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 c1, c2, 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. 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—6to 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 Sect. 2, 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 Sect. 3, 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 Sect. 4.

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:

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. [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.

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.

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 Pfwas 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 Pfwas full, the failed individuals in Pfwere updated based on the‘‘first in,first out’’rule.Pfmaintained 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., Fi～N（uF，0.1）, which can yield better results compared with the Cauchy sampling method, and uFwas updated [28] at the end of each generation as follows:

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

where SFis a set for all failed mutation factors in each generation, and the uFsaved in SFis used to provide the prior historical information for the next uFupdate. SFis 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 uFwith generations from 5 to 20 [28];uF=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:

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 CRiis generated,as shown in Eq. (9),and updated using Eq. (10) [28].

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

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:

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 follows:

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 106as the iteration termination criterion to replace the maximum iteration number; this reduced the operating time and improved the quality of the final optimal solution.

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:

Fig. 3 (Color online) Relationship between QS and fitness of ‘‘for BNCT’ spectrum

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.

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 106. The average termination generation for 20 runs of the ‘‘Cf-252’’ spectrum was the 1157th generation; ‘‘Cf-252 Mod’’ spectrum,the 544th;‘‘for BNCT’’spectrum,the 451th;and‘‘reactor’’spectrum, the 472th. 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. 5 (Color online) QS of solutions with average over multiple independent runs in SDENUA

Fig. 6 (Color online) Termination generation of final optimal solution of 20 runs when fitness of solutions reached 106

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.

The oscillations in the unfolded spectra are an inherent feature of the numerical solution to the Fredholm equation[19, 20, 31], 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.

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 (Cujand Ruij) were generated using Eqs.(14)and(15),respectively,and the chi-square was set to 10-6to compare with the fitness upper bound in the SDENUA. Figure 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—9to 10-7MeV. As shown in Fig. 2, the response functions of energies from 10—9to 10—7MeV overlapped significantly,thereby weakening the unfolding power of the response functions in this energy range.Owing to the integral fitness function,as shown in Fig.8c,d,although the final optimal solutions fluctuated in a large range around the reference spectra, the fitness of the final solutions still reached 106.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.

Fig.7 (Color online)QS of final solutions from SDENUA smoothed using Gaussian and without smoothing

Fig.8 (Color online)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,

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.

Fig.9 (Color online)Uncertainties of unfolded results.Uncertainties for calculated neutron counts were less than 5% (calculations are presented in the ‘‘Appendix’’)

Although the fitness value of the final solution was limited to 106, 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 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 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 106.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,106as 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 106rapidly 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%.

Appendix

In this study, the type A standard uncertainty [33] was used and calculated as follows:where ΔUAis the uncertainty of type A with the same units as variable x,xithe calculated value of variable x based on the spectra uncertainty estimation (xiis the calculated neutron fluence of the ith run; in the neutron count uncertainty estimation,xiis the calculated neutron count of the ith run),the mean of variable x, and n the number of independent runs. The uncertainty of type A in percentage is estimated as follows:

where xexacis the exact value of variable x.