Optimization of load-follow operations of a 1300MW pressurized water reactor using evolutionnary algorithms

Because of the increase of intermittent renewable energies, load follow operations for French PWR will be crucial in the years to come. The goal of this study is to make realistic changes to the plant operations in order to optimize load-following without disrupting the plant. 6 discrete parameters were considered among the overlaps, speeds and maneuvering bands of the control rods. A simulator oriented model of the pressurized reactor based on APOLLO3 R © is used. It includes 3D neutronics calculations with point kinetics and a 0D model of the secondary system. The operating mode (G mode) was modeled, to account for the human operator of the power plant. Two objectives were considered so as to both minimize the volume of effluents generated during the transient, and maximize the core axial stability. The reaction of the power plant to a load follow transient varies greatly during the operating cycle, because of fuel depletion effects. Therefore, 4 burnup points are considered and the objectives are computed for each point and then reduced to two “whole cycle” objectives. A biobjective massively parallel asynchronous master worker evolutionary algorithm based on AMS-MOEA/D was implemented. It is a highly exploratory algorithm, suitable for black-box problems with an important computing time. The analysis of the performances of this algorithm shows that it is able to find a diversified Pareto front, with solutions that greatly improve the load follow operations for all burnup points compared to the standard control rod parameters.


INTRODUCTION
Through the combined effect of the increase of renewable energies and the shutting down of highly carbonated power plants, France's energetic transition faces an issue of balance between consumption and production [1]. French agency RTE describes an electricity market without margins, and highly dependent on nuclear power plant (NPP) performances, because of the prominent place of nuclear power in the production mix. With a future mix made of 50% of nuclear energy and a high portion of renewable energy (at least 30%), nuclear power faces a manageability challenge. NPP are designed so that they are self-regulating in the case of small power variations. However, in load-following transients, shifts in power can be important and therefore temperature can be quite far from its nominal value. Therefore, the management of control rods and soluble boron is crucial to mitigate temperature variations. Moreover, axial instabilities can develop because of the evolution of poisons such as xenon, and need to be mitigated by the operator. This work aims at optimizing the response of a NPP during an electrical power transient according to some values of interest by changing the control rod configuration. The values of interest considered are the volume of effluents generated by the transient, directly related to the cost of the operation, and the axial power imbalance, related to the safety of the operation, and the ability for the reactor to undergo a new transient without having to wait for the end of a xenon oscillation. Because of fuel depletion effects, one rod configuration will not maintain a constant performance during the operating cycle. It is therefore necessary to find configurations that will be optimized for the complete cycle. A simulator oriented model of a French 1300MW PWR was developped and is used in this study [2], especially targeting load-following transients. The core calculation are performed thanks to APOLLO3 R [3], with a model including 3D neutronics, multi1D thermal hydraulics and fuel thermal effects, and a 0D kinetic model. This model is coupled to a 0D model of the secondary system. The NPP is controlled via a simulated operator based on the current G mode. The computation time for one burnup point is around 10 minutes.
The search space for the rod configuration has a size of 10 12 , therefore a strongly exploratory algorithm is needed for the optimization process. A massively parallel MOEA/D [4] biobjective algorithm was implemented, close to the one developped in [5], and adapted for the mutliple burnup problem. A fitness landscape analysis [6] was first conducted on the search space by computing a random walk on the parameters in order to tune the individual-burnup algorithm parameters as was suggested in [2]. The individual-burnup optimization of load follow operations at 4 burnup points was successfully conducted, using the results from this random walk.
In this paper, a methodology for the multiple-burnup (MBu) optimization of load follow operations is presented. Its goal is to minimize the cost of a load follow transient by reducing the volume of effluents, and to ensure maximum safety by minimizing a function representative of the integrated axial power imbalance during the operation at 4 burnup points simultaneously. First, the case study and simulator are presented, then, a massively parallel evolutionnary algorithm based on the MOEA/D algorithm is developped. Finally, the results of the method are compared to the individual-burnup (IBu) method, and the solutions are discussed.

Description of the Reactor
The reactor studied in this work is a 4-loop French 1300MW PWR. Its maneuvering potential is important, and higher than that of the French 900MW PWR [7]. The core is 4m high, cylindrical and contains 193 assemblies, 120 made of Uranium oxyde (UO2) and 73 made of Uranium plus Gadolinium oxides. The reactor fuel is managed within the GEMMES mode, allowing an extended irradiation cycle, by the use of natural Uranium enriched to 4% in UO2 assemblies and a fuel reloading by third.
The reactor is conducted with the "G mode". Two kinds of control rods are used in this mode,

Transient Simulator
The simulator used for this project is a multi-physics PWR model encompassing both the primary and secondary systems [2] [8]. The computations are performed thanks to APOLLO3 R [3], the multi purpose deterministic transport code developped at CEA. The neutronics model is based on a static 3D flux calculations (two-group diffusion equation) for the power distribution calculation, and a point kinetics model (linearization of the Nordheim's formula) for the core power calculation. This model takes into account Doppler and moderator feedback effects through the use of a simplified thermal-hydraulics model computed in APOLLO3 R . The evolution of actinides and several poisons, including Xenon, Iodine and Samarium is computed at each time step. This neutronics model is coupled to a 0D model of the steam generator, so that the simulator can take into account the imbalance between the extracted steam flow and the produced steam flow during loadfollow operations. A simplified model of the G mode is used in order to reproduce rod movements and boron concentration evolutions [2].

. Optimization Problem
The purpose of this simulator is to compute the objectives of the bi-objective optimization study of load-follow operations. The considered operation is a classic day/night load-following scenario (see Fig. 2). The power is set to 30%NP (Nominal Power) for 6 hours and then goes back to 100% NP. The slope is the maximum authorized slope of 5% NP/min. The optimization is biobjective, with two opposite goals: assuring the axial stability of the transient, by minimizing an axial-offset integral and reducing the cost of the transient by minimizing the volume of boron effluents. Safety is also considered via constraints related to the anti-reactivity margin. For this optimization problem, the free parameters are the maximum speeds of the PSR rods (2 discrete parameters S1 and S2 * ), the overlaps between the PSR rods (3 discrete parameters O1, O2 and O3) and the maneuvering band of the TRR (1 discrete parameter M B). Their ranges are limited by technological and safety constraints, and listed in Tab. 1. 4 burnup points are considered throughout the cycle: 0% of the cycle length, 20%, 50% and 80%.

ALGORITHM
Several approaches are available to search for globally interesting solutions. First, one simple approach is to optimize the operation at one burnup point with the IBu method, and then compute the results at the other burnup points, in order to see if they maintain their performance. This is motivated by the fact that these solutions have interesting properties, and thus a higher chance than the general population to perform well at other burnups. This technique was used by Dr. Muniglia in [8] for this specific problem. However, there is no formal proof that the results of this method will be optimal. Therefore, an approach that tackles the 4 burnup points simultaneously is necessary.
It is possible to consider each objective at each burnup as an objective, transforming the problem into a 8-objective optimization problem. Algorithms exist in the litterature to solve this type of optimization problems [9]. However, they are difficult to implement. Moreover, this would lead to a 8-dimension decision surface, making it hard for a plant operator to choose one configuration. In order to facilitate both the decision making and the optimization problem, we choose to reduce the dimension of this problem to a biobjective problem, by computing two aggregated objectives from the calculations at each burnup. With v(bu) and i(bu) the volume of effluents and imbalance integral at one burnup point bu, the multiple-burnup (MBu) objectives are as follows: The MBu volume of effluents is computed as the weighted mean of the 4 IBu volumes. The volume of reference for the weights is the volume computed for the base parameters. This is motivated by the fact that the volume of effluents is higher towards the end of the cycle, because the boron concentration is more important at the beginning of the cycle. This choice of weights aims at giving the same importance to the 4 burnup points. The MBu imbalance is computed as the maximum of the IBu imbalances. This is motivated by the fact that imbalance is first a safety concern, and therefore the worst case scenario should be considered.
The two MBu objectives are then used in the MOEA/D bi-objective algorithm. The MOEA/D algorithm is based on the decomposition of the bi-objective problem into λ mono-objective problems, cutting the objective space along λ directions. For one direction i, the fitness is computed as: with w i = (cos (iπ/(2λ)) , sin (iπ/(2λ))) the weights associated with direction i, and z * the objective point. Each direction searches for an optimal solution for this scalarized function with an (1+1)-evolutionary algorithm which iteratively mutates and updates the best known solution. In order to converge to a diverse Pareto set, directions can exchange solutions. In our case, as suggested by Dr. Muniglia in [8], we use λ = 200 directions, with no limitations to the communications between directions. The mutation operator was tuned using fitness landscape analysis. The algorithm is implemented as an asynchronous master-worker system in order to face the large variance of computation time between fitness evaluations [10]. The objectives are normalized using standard normalized variable, based on the data acquired from the fitness landscape analysis [2].

RESULTS
The   most of the solutions obtained with the MBu optimization method dominate the solutions found by the IBu optimizations, confirming the interest of the method. In particular, these results confirm that optimizing only at 0% of the cycle length is not enough to ensure that the plant will perform well for the rest of the cycle. Moreover, the MBu managed to dominate the other method with a very reduced computing time. Fig. 4 shows the results of the MBu method for each burnup points. Overall the results of the MBu method present good performance on the individual burnup points, with points close to the reference Pareto front, except at 0% of the cycle length. This gap between the two methods for this burnup point is likely due to the fact that the variance of the volume of effluents at the beginning of the cycle is lower than at the other burnup points. Therefore the MBu method has less interest in reducing the volume of effluents at this point since it will not impact the MBu objective as much.
Finally, Fig. 3b shows the solutions in the parameter space. The algorithm's result is a diverse Pareto set, therefore the solutions are diverse. However, some general conclusions can be drawn. First, the maneuvering band of the TRR should be extended, giving more room for the operator to compensate the xenon transient without using a lot of soluble boron. This, however, intuitively worsens the power imbalance, since the TRR is made of black rods operating in the upper part of the core, pushing the power towards the bottom of the core. By reducing the overlap between the second and third groups of the PSR, and by reducing the speed of the first group, it appears to be possible to counter that effect and achieve a gain on both the power imbalance and the volume of effluents. These compensation effects between the PSR and TRR are the results of shadow effects between the different rods, and are not intuitive, therefore confirming the interest of a retroengineering method based on optimization. Overall this indicates that the TRR should be given more importance in the mitigation of temperature differences, while making sure that they do not create important axial imbalance. It can be noted that the "T Mode" on the European Pressurized Reactor, designed to optimize load-follow operations, does have a more lenient way of managing the control rods [11].

CONCLUSIONS
In order to optimize the availability of nuclear power plants for complicated load follow transients with important slopes and amplitude, optimizing the rod configuration is a strong lever, allowing the operator to reduce both the axial imbalance resulting from the transient, and the volume of effluents, hence the cost of the operation. The method proposed here underlines the fact that optimizing the operation at one point of the cycle is not enough to guarantee good behaviour on the whole cycle. A multiple-burnup method was proposed, giving better results than individualburnup optimization, with a reduced computing time.
Solutions were found that perform well and even close to optimality at each considered burnup points. It was found that these solutions are not intuitive, therefore confirming the interest of the optimization process. One key result of the method is that the maneuvering band of the Temperature Regulation Rods should be extended, and the overlaps between the Power Shimming Rods adapted. Furthermore, the method can be adapted by changing the multiple-burnup objective definition, in order to better match the needs of a NPP operator.
Though the algorithm already gives satisfying results with a reasonable computing time, further work will tackle the acceleration of its convergence with metaheuristics in order to improve the algorithm's performance so that it can be used on a variety of load-follow transients with a reasonable computing cost.