Multiobjective Core Reloading Pattern Optimization of PARR-1 Using Modified Genetic Algorithm Coupled with Monte Carlo Methods

Center for Mathematical Sciences (CMS), Pakistan Institute of Engineering & Applied Sciences, Nilore 45650, Islamabad, Pakistan Department of Nuclear Engineering, Pakistan Institute of Engineering & Applied Sciences, Nilore 45650, Islamabad, Pakistan College of Engineering, King Saud University, PO-Box 800, Riyadh 11421, Saudi Arabia National Tokamak Fusion Program, Nilore 45650, Islamabad, Pakistan


Introduction
Nuclear fuel management is one of the most significant considerations when managing smaller research reactors. In these applications, the reactors must secure a wide variety of safety restraints while also providing high thermal neutron flux for external experiments. Core optimization can be performed at the assembly level by optimizing fuel composition or at the core level by optimizing the assembly loading pattern. Overall, many aspects of core operation are very sensitive to loading patterns such as power density, burnup, reactivity distribution, and accumulation of fissile materials. Naive calculations for core loading optimization techniques can be prohibitively expensive due to the computational cost of simulating each arrangement of the reactor core. is is especially true because a design may need to simultaneously optimize multiple parameters such as peaking power, radial leakage, and end-of-cycle core reactivity [1]. is paper presents an intelligent technique using modified genetic algorithms to perform these tasks with reasonable computational cost.
Extensive studies have been done to optimize the in-core fuel loading pattern during the last few decades, resulting in numerous perspective processes. e loading pattern must ensure the safety and efficiency of a research reactor. Several different techniques were used in the past for core loading optimization; some of them were simulated annealing method [2][3][4], genetic algorithms (GAs) [5][6][7][8], simple tabular search [9], particle swarm optimization [10], and artificial neutral network [11][12][13][14]. e GA was first introduced by Holland in the 1970s as one of the stochastic optimization techniques inspired by the mechanism of the natural evolution [15]. It has demonstrated to overcome traditional techniques in many complex problems. By using a GA, an in-core fuel management strategy has been developed in the past to optimize more than one parameter [2]. e GA differs from most of the optimization techniques by searching for one group of solutions to another, rather than focusing on one solution to another.
is makes it suited uniquely to most of the multiobjective optimization problems. An improvement of new coding procedure for GA is capable of searching automatically for the optimal fuel loading patterns most suitable for the research reactor [15]. e standard bit-based genetic operators in GA are used to optimize the arrangement of assemblies, burnable absorber, and burnt assembly orientations in the pressurized water reactor (PWR) [12]. e application of GA to reloading pattern optimization problems of a PWR allows pattern from different distributions to combine and reproduce. An optimized pattern in GA is better than that obtained from linear programming [14]. e multiobjective GA (MOGA) used three integer-based arrays to perform the multiobjective optimization of fuel loading patterns of PWR. e MOGA makes it possible to identify all the competing objective surfaces in a single optimization run [5]. In comparison of GA to simulated annealing, the GA optimization performed better in the initial global search, but simulated annealing was better for the local search. GAs seem to perform better in solving LP optimization problems [14].
In nuclear design calculation, the effective multiplication factor k eff and the neutron flux distributions are considered the most foundational evaluation quantities. It is used broadly to derive many important formulations in nuclear design calculation such as control rod worth and the reactivity coefficient. Optimizing the placement of fuel assemblies in the core is the most effective way to increase the utilization of a research reactor. erefore, one needs to maximize k eff and neutron flux densities inside the reactor core [16][17][18][19][20]. In the present study, loading patterns for PARR-1 have been proposed by maximizing the effective multiplication factor and neutron thermal flux while keeping low power peak. e main objectives of the fuel management optimization problems considered in the present work are listed below: (i) Maximize the k eff (ii) Maximize thermal neutron flux in the flux trap e study follows the two constraints: (i) Keep the local power peaking factor lower than a predetermined value to maintain fuel integrity.
(ii) e positions of control fuel elements (CFEs) in the core are fixed, but they can swap their position among the CFEs.

SuperMC Computer Code.
Practical computer simulations of a nuclear system like reactors use either deterministic or Monte Carlo codes. Due to continuous energy cross-sections and flexible geometrical capabilities, the Monte Carlo (MC) method is distinctive to simulate complicated nuclear systems and is envisioned as a routine method for nuclear design and analysis in the future. e SuperMC computer code, developed by the FDS team in China, is a general-purpose CAD-based MC-based radiation transport code. It can simulate neutrons, photons, and their coupled behavior. It performs transport calculation, geometry and physics modeling, and visualization of results and geometry. It has been developed and verified by using a series of benchmarking cases such as the fusion reactor ITER model and the fast reactor BN-600 model. SuperMC is still in its evolution process toward a general and routine tool for the simulation of nuclear systems [21]. SuperMC computer program exploits the hybrid MCdeterministic method concept and advanced computer technologies. Being general-purpose radiation transport program, SuperMC is designed for high-fidelity simulation of nuclear-system problems such as reactor analysis, radiation shielding and dosimetry, medical physics, and radiation detection. SuperMC code can be applied to transport calculation of various types of particles, depletion and activation calculation including isotope burnup, material activation and shutdown dose, and multi-physics coupling calculation including thermo-hydraulics, fuel performance, and structural mechanics. e powerful bidirectional automatic conversion between general CAD models and calculation models can be easily performed. Results and the process of simulation can be visualized with a dynamical 3D data set and geometry model.
Advanced cloud computing framework makes the simulation, which is extremely intensive in computation and data storage, more attractive just as a network service. e modular design and generic interface promotes its flexible manipulation and coupling of external solvers. e architecture of SuperMC is shown in Figure 1.

General Description about the Strategy Developed.
e optimization strategy presented in this paper is aimed at producing the best possible core loading pattern for any core by modeling of the core using SuperMC. It is based on the probabilistic approach of genetic algorithm. e detail of the application of the genetic algorithms will be given later in the description of the code. e variables normally available in selecting a fuel loading plan are fuel enrichment of fissile content of the fuel, the number of fresh fuel assemblies to be loaded, the arrangement of the fresh and partially spent fuel assemblies in the reactor, and the techniques used to control the excess reactivity of the reactor during the cycle. e code has provision for optimization based on three different factors, i.e., it optimizes the core on the basis of maximum reactivity, which in turn represents longer cycle length; maximum thermal or fast flux in a fixed element, e.g., flux trap or beam tubes for material test reactors; or maximum flux in a particular region of the core. e optimization of the core is also subject to a constraint on the power peaking factor. e value of power peaking factor must not exceed the designed maximum limit for the reactor. is has been incorporated by reducing the fitness of the loading pattern that violates the constraint so much that chance of its propagation becomes almost zero. e user also has the liberty to select the population size, the total number of generations, crossover fraction, mutation fraction, and the fraction of population of loading patterns that is to be promoted into the next population as being elite.

Description of the PARR-1 Core.
e reactor core considered in this study is Pakistan Research Reactor-1 (PARR-1), the swimming pool-type material test research reactor (MTR). PARR-1 has a full power of a 10 MW, utilizing the low-enriched uranium fuel (LEU) of uranium silicide (U 3 Si 2 -Al) fuel contained less than 20% of enriched 235 U. e reactor core consists of identical standard fuel elements (SFEs), control fuel elements (CFEs), and an irradiation water box (WB) at the center of the core created by replacing a fuel element. e core immersed in the demineralized light water acts as a coolant, moderator, and reflector. e core is reflected on two sides left and right by the light water, while behind a gamma lead shield of thickness 12.7 cm, light water, and graphite column reflect the third side. Moreover, both graphite and light water reflect the fourth side [22]. e core filled with assemblies with FA types is shown in Table 1. Core configuration modeled by SuperMC and its operating conditions are given in Figure 2 and Table 2, respectively.

Description of PARR-1 Fuel Assembly.
e full core consists of 6 × 5 arrays, containing 29 fuel elements and a flux trap in the center and assembled with partial graphite reflector and one thermal column, all of which are shown in Figure 2. e standard of the fuel element is modeled in three regions. Region 1 is home to 6.275 cm × 8.10 cm uranium metal. Regions 2 and 3 are thin claddings of 0.718 cm × 8.10 cm aluminum in Figure 3(a). Figure 3(b) represents the control fuel element, which is divided into nine regions. Region 1 located in the central part housing an absorber rod with dimensions of 6.275 cm × 3.370 cm. Twofuel regions 2 and 3 are above and below this absorber region with a dimension of 0.718 cm × 2.1965 cm, respectively. Six regions on both sides represent the side plate with dimensions of 0.718 cm × 2.5335 cm, 0.718 cm × 3.370 cm, and 0.718 cm × 3.1965 cm, respectively. e material composition of the control rod is 80% Ag, 15% In, and 5% Cd. All these five control rods are used as safety rods, one of which with the lowest reactivity worth is used as a power regulating rod [22]. e standard fuel element and control fuel element modeled in SuperMC are shown in Figure 4.

Core Neutronic Calculations Using SuperMC.
e core neutronic calculations for PARR-1 to evaluate the objective functions, i.e., k eff and thermal fluxes were preceded using SuperMC. e simulations for the PARR-1 core neutronic calculations were run for 100,000 neutron histories for a total of 120 cycles (20 inactive cycles, 100 active cycles). e 100 independent simulations were run to get the objective function values corresponding to the 100 randomly generated different loading patterns as an initial population of MGA.

Modified Genetic Algorithm Characteristics.
To determine, which chromosomes will be chosen as the basis of the   Science and Technology of Nuclear Installations 3 subsequent generation, the reproduction is considered in the simple form of conventional GA. e best chromosome from the last population might be lost from generating populations from only two parents. e selection of crossover or mutation or both operations to generate new chromosome may not reach the best fitness of the objective function value. ereby, it can be said that the best fitness found through the simple genetic algorithms cracked from the new population may be inferior to the old generations. e objective of the proposed MGA is to avoid this disadvantage. e structure of the MGA is quite similar to the simple GA. However, the MGA has been distinguished from the simple GA in that the new offspring are generated after both the crossover and mutation operations have been applied. us the deterioration problem will not take place on the grounds that the best fitness from the new generated offspring will be the fittest or at least the similar with the previous fitness.
ere are many operators, parameters, and other settings involved in the genetic algorithms. ese ingredients for a GA can be incorporated differently in various problems. e appropriate selection of these ingredients in the evolutionary process plays a vital role in convergence; otherwise, the process will make the GA susceptible to premature convergence. erefore, the better choice of a specific crossover depending on the nature of the problem can improve the performance of the genetic algorithm. e genetic parameters like crossover and mutation probability, numbers of individuals and generations should be determined depending on the regarded problem. e population size and the number of generations are particularly considered to be the important parameters of GA. In case of too low number of chromosomes, the crossover operator will have a less number of opportunities to perform crossover. en, it will explore only a small part of search space. On contrary, if there is too large number of chromosomes, the performance of the GA will be slowed down. Several runs of MGA have been    Table 3.
Optimization of PARR-1 using MGA and SuperMC is the novelty in this paper. e core of PARR-1 has never modeled before using SuperMC code and never been optimized using MGA in the past. e name-modified GA has been used because it is different in some aspects from a traditional GA. In a conventional GA, reproduction in subsequent generation is considered to choose the chromosomes for next generations. e best chromosomes from  Science and Technology of Nuclear Installations the last population may get missed by generating population only from two parents. e selection of crossover or mutation or both operations to generate new chromosomes may not converge to the best fitness value of the evaluation function. is disadvantage of conventional GA can be avoided by the application of proposed modified GA. In the modified GA, the new descendants are generated by the application of both crossover and mutation operations. us the missing of the chromosome having better fitness value than the previous can be avoided in this way or at least the offspring will have the similar fitness than the subsequent generation.
e flow scheme of the designed strategy is given as follows: (1) e total number of assemblies; different types of assemblies; their indices; the optimization method; and limits on the constraints, the optimization function, and the assemblies are saved along with the parameters required for genetic algorithms from the file "CORES.inp." (2) e Python code was developed to generate the SuperMC readable input file with the changing core loading pattern for each simulation.  is process continues until the maximum allowed number of iterations has been reached or the given tolerance level is met. e flow chart in Figure 5 depicts the whole process.

Numerical Results and Discussion
e optimization of loading pattern for the PARR-1 core was performed subject to a constraint on the power peaking factor. e fitness function is defined for obtaining the best configuration of the reactor to obtain the maximum k eff under safety considerations of the maximum allowable power peaking factor of 2.089 [16]. is has been incorporated by reducing the fitness of the loading pattern that violates the constraint so much that chances of its propagation become almost zero. Initially, 100 randomly generated loading patterns were simulated as an initial population. e total number of generations, crossover fraction, mutation fraction, and the fraction of population of loading patterns that was to be promoted into the next population as being elite were set for the calculation of the best loading pattern and are presented in detail in Table 4. e total number of assemblies, types of assemblies, their indices, limit on the power peaking constraint, the optimization function, and the assemblies were input along with the parameters required for GA. All the possible locations of each assembly type were noted. e indices assigned to each assembly were saved in a one-dimensional array that served as the chromosome for the GA. e initial population is generated of such different chromosomes containing these indices at different locations. For each chromosome, an input file is generated using the Python script, which is readable for SuperMC in which the position of each assembly depends on the position of its corresponding index in the chromosome. e objective function values for k eff and thermal neutron fluxes in the flux trap were obtained using SuperMC. e values of the optimization functions were saved. e power peaking factor was calculated. e chromosomes resulting in a peaking factor greater than the limiting value in the core are placed on the least priority. Parents for crossover are selected using the function to assign high probability to the chromosomes with high fitness. A new population was generated using GAs and the defined fractions of different ways of generation like mutation, crossover, and elite count. A total of 100 generations were simulated in a similar way to get the optimal solution.   Table 3 shows the sensitivity of initial population.
100 numbers of independent runs were made, and it is observed that 100 numbers of individuals give optimized results. It is also observed that any increase in population size would increase computational cost without increasing   the accuracy of results and would be unnecessary. Based on this sensitivity analysis, initial population size of 100 is selected for further analysis in this research.

Stability Analysis of MGA.
e test of the algorithm stability is one of the most important tests that should be conducted. e objective to check the stability of the algorithm is the issue that the algorithm gives the same answers, and the uniqueness of the optimized answer should be tested. To this end, GA was run 100 times independently and it was observed that the results are in 95% confidence interval. e results of the first 10 repeated runs of the algorithm along with the mean value and the variance are shown in Table 5 and Figure 1.
Results shown in Table 6 and Figure 6 apparently represent a small difference among results of repeated simulations. e 0.0000003 variance in k eff and the 0.00002 variance in maximum power peaking factor are very low, which shows high stability of the modified genetic algorithm in different runs for 100 repetitions. e convergence of the objective function values corresponding to the best loading patterns in each generation    with the maximum k eff value and maximum thermal neutron flux is illustrated in Figures 7 and 8, respectively. e maximum power peaking factors for the cores during the convergence are also shown in Figures 7 and 8 Tables 6 and 7 show the groups of generations having the same optimum value of k eff , thermal fluxes, and the maximum value of power peaking factors of th core, respectively. e best optimum loading patterns for the group of generations having the same value of k eff and thermal fluxes are shown in Figures 9 and 10, respectively. e colormap shows the normalized radial power distribution in the core.
Summary of the computational time to get the best optimal loading patterns is given in Table 8. e parallel simulations of the Monte Carlo code for the calculation of objective function values were run on 50 CPUs. It takes approximately 122 core-hours (CH) to complete 100 independent simulations on 50 CPUs. To find the best optimal loading pattern for PARR-1 core, the GAs take around ∼0.1 core-hours (CH) on the master node with system specifications of Intel (R), Core (TM) i7-7700 CPU @ 3.60 GHz 3.60 GHz. erefore, the total computational time to get optimal solution was approximately 122.1 CH.

Conclusions
is study has demonstrated a powerful multiobjective optimization method for the nuclear fuel loading optimization problem of the MTR-type research reactor PARR-1.
e SuperMC code has been successfully applied for the neutronic calculations of PARR-1 core and combined with the modified GA to evaluate the objective function values. e results of optimization illustrate the efficiency of the MGA search for the optimal fuel loading pattern to obtain its maximum effective multiplication factor and the maximum thermal neutron fluxes in the flux trap region while keeping the power peaking factor under design limits. It is worth mentioning the results that the optimum value of k eff is 1.038659 and the optimal thermal neutron flux in the trap is 7.93501 × 10 13 n′s/cm 2 /sec. Moreover, based on the present research it can be concluded that MGA coupled with SuperMC can be efficiently used for the loading pattern optimization of PARR-1. Also, the best optimal loading pattern can be obtained in around ∼122.1 core-hours on the previously mentioned computer system details.

Data Availability
Due to confidential matters, data might not be shared.

Additional Points
(i) Optimal core reloading patterns were generated for Pakistan Research Reactor-1. (ii) SuperMC has been used for the core neutronic calculations. (iii) A powerful multiobjective optimization method GA was coupled with Monte Carlo method. (iv) k eff and thermal neutron flux were considered as an objective function for optimization. (v) e local power peaking factor was kept lower than a predetermined value. (vi) e best reloading patterns were recommended by using the implemented strategy.

Conflicts of Interest
ere are no conflicts of interest.