Baptiste Gasse

Sébastien Verel
email: verel@univ-littoral.fr

Jean-Michel Do
email: jean-michel.do@cea.fr

MOEA/D with Adaptive Mutation Operator Based on Walsh Decomposition : Application to Nuclear Reactor Control Optimization

Keywords: Applied Computing Methodologies, Bi-objective Optimization, Surrogate model/fitness approximation, Nuclear Reactor Physics

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION 1.Context and Motivations

In France, the strategy adopted in terms of energy transition is leading to a reflection on a change of the national electricity production mix to achieve carbon neutrality by 2050. In 2022, 62.7% of electricity production was provided by nuclear reactors in metropolitan France, 11.1% by hydroelectricity sources, 12.7% by wind and solar power, and 13.5% by fossil-fired power stations, the vast majority of which are gas-fired. While the energy produced by the fleet of nuclear reactors is controllable, i.e. capable of adapting their power output to ensure the supplydemand balance of the electricity network, wind and solar energy sources are intermittent since weather conditions affect their load factor. This technological difference coupled with the fact that the merit order curve of the French electricity market gives consumption priority to intermittent resources raises a growing need for manoeuvrability for current nuclear reactors.

Indeed, research has shown that a high penetration of intermittent energies in the electricity production mix decorrelates the daily nuclear production profile from real electricity demand, implies reduced power variation times and, depending on the number of reactors in operation in the coming years, decreases the nuclear fleet's load factor [START_REF] Denholm | Grid flexibility and storage required to achieve very high penetration of variable renewable electricity[END_REF][START_REF] Cany | Nuclear contribution to the penetration of variable renewable energy sources in a french decarbonised power mix[END_REF]. Numerous works have studied the issue from the point of view of economic constraints, leading to unit commitment type optimization calculations in electricity mix scenarios with a high penetration of intermittent energies without relying on fossil plants [START_REF] Xu | Unit commitment model considering nuclear power plant load following[END_REF][START_REF] Jenkins | The benefits of nuclear flexibility in power system operations with renewable energy[END_REF][START_REF] Ju | Unit commitment accommodating large scale green power[END_REF], but also with technical feasibility criteria since power variations in a pressurized water reactor (PWR) induce operating constraints linked to safety rules. This article introduces a work that is in line with this second approach.

The optimization problem at stake in this study requires the use of a black-box function with a high computation time, as around 45 minutes are needed for the calculation of a single point in the objective space. In a context of limited computing resources, this time constraint implies that the number of evaluations can't be too high even though the main goal is to compare optimization methods with each other while guaranteeing sufficient statistics to interpret the results. The use of best direction search assistance in the decomposition of the search space for the MOEA/D algorithm using a surrogate model of the fitness function has led to a slight gain in performance [START_REF] Drouet | Surrogateassisted asynchronous multiobjective algorithm for nuclear power plant operations[END_REF], but the use of an adaptive mutation operator for a reactor controlling simulator hasn't been tested yet. In this study, we propose a new adaptive method based on a Walsh surrogate model of the fitness function to adapt regularly the mutation operator of MOEA/D.

The paper is organised as follows. Firstly section 2 introduces related work in the fields of MOEA/D adaptive methods, use of Walsh surrogate models for optimization techniques and nuclear reactor optimizations for load following control mode. Section 3 is reserved both for a review of nuclear reactor physics notions useful for the understanding of PWR control, and the definition of the optimization problem addressed in this study. Section 4 describes the methodology behind the adaptive bi-objective evolutionary optimization algorithm using Walsh surrogate model to update the mutation operator regularly. Finally, section 5 describes the experimental setup and some discussions of the results obtained with our method.

Related Works

Adaptive evolutionary algorithms : in [START_REF] Shim | A hybrid adaptive evolutionary algorithm in the dominationbased and decomposition-based frameworks of multiobjective optimization[END_REF] a general optimization study using three evolutionary algorithms of different natures has been conducted. Authors found that using the ratio of promising solutions obtained by each optimizer at a given generation to select the proportion of distributed solutions for each algorithm at the next generation results in a significant performance improvement compared with non-adaptive techniques. Adaptive weight adjustment techniques for MOEA/D has been conducted in [START_REF] Qi | MOEA/D with Adaptive Weight Adjustment[END_REF] in the case of complex Pareto front to prevent the penalization of some regions in the objective space during the optimization. More recently, optimization methods involving an adaptation of the mutation operator in an evolutionary algorithm has been carried out in [START_REF] Cruz-Salinas | Selfadaptation of genetic operators through genetic programming techniques[END_REF][START_REF] Prieto | Hybrid adaptive evolutionary algorithm for multi-objective optimization[END_REF], where decision trees and decomposition-based approaches helps to adapt dynamically the mutation rates of genetic operators. This multi-objective hybrid adaptive algorithm has shown better performance on benchmark functions in its ability to find a well-covered and well-distributed set of points on the Pareto Front.

Black-box optimization : in [START_REF] Dadkhahi | Fourier representations for black-box optimization over categorical variables[END_REF], two novel Fourier representations as surrogate model for problems involving the use of black-box functions over categorical variables has shown good performance when combined with a variety of search algorithms. One of the Fourier representation developed in this study is made of an abridged one-hot encoded Boolean expansion recalling previous work published in [START_REF] Verel | A surrogate model based on walsh decomposition for pseudo-boolean functions[END_REF][START_REF] Leprêtre | Walsh functions as surrogate model for pseudoboolean optimization problems[END_REF][START_REF] Derbel | Walsh-based surrogate-assisted multi-objective combinatorial optimization: A finegrained analysis for pseudo-boolean functions[END_REF] which developed Walsh polynomials as surrogate models for pseudo-Boolean optimization problems. These studies have shown Walsh surrogate models are relevant for benchmark problems and should be considered to improve the performance of optimization algorithm for combinatorial problems.

Nuclear reactor control optimization : previous works [START_REF] Muniglia | A fitness landscape view on the tuning of an asynchronous master-worker ea for nuclear reactor design[END_REF][START_REF] Drouet | Surrogateassisted asynchronous multiobjective algorithm for nuclear power plant operations[END_REF], on optimising the control of nuclear reactors have shown the existence of optimized sets of technical parameters in order to minimize the axial destabilization of the power profile within the reactor core (see sect. 3) as well as the use of liquid effluents in load following mode. These optimizations were carried out for a fixed transient and by exploring different parameters for the movement of the control rods which are inserted into the core during the power reduction. In our study, we propose to jointly optimize rod parameters with the best way of distributing a fleet transient involving two nuclear reactors.

NUCLEAR REACTOR PHYSICS

Description of a French PWR reactor

This optimization study is conducted on a fleet of two 1300 MW nuclear reactors typical of the units involved in load following of the French power grid. Apart from the heat production medium, a nuclear reactor does not differ much from a conventional thermal power plant in terms of electricity production.

There are three circuits : first, the primary circuit, where water passes through the reactor's core from bottom to top and is heated by the fission reactions of uranium 235 nuclei. Once it reaches the top of the core, the primary water is redirected to a steam generator where its heat is exchanged with the water of the secondary circuit, which vaporizes. The steam thus produced is expanded in the turbines, where thermal energy is converted into mechanical energy. The rotation of the shaft of the successive turbines makes an alternator to rotate, which converts kinetic energy into electricity. The role of the last circuit, called the cooling circuit, is to evacuate the residual heat which could not be converted into mechanical energy by virtue of the second principle of thermodynamics. The core is roughly comparable to a cylinder 4 meters high and 3 meters in diameter. It is composed of 193 fuel assemblies and 62 reflector assemblies (see fig. 1, left). The fuel assemblies all have the same geometric pattern : 264 fuel pins about 4 meters long that contain uranium oxide pellets (UOX) over their entire height, and 25 guide tubes (see fig. 1, right). These tubes enable neutron-absorbing rods to slide vertically in order to regulate the neutron population, and thus the fission chain reaction. These rods are placed at strategic points in the core for the sake of a flatten radial flux distribution during power variations. The reflector assemblies, shaded in gray on the figure 1 (left) limit the leakage of neutrons out of the core as well as the wear of the vessel which ages with irradiation.

Rods

Concerning the rods, the rate of captured neutrons increases as they inserts into the core. Transiently, the number of neutrons decreases in the core and consequently the fission rate, which implies a decrease in the power produced. The opposite occurs when rods are extracted from the core. Rods are divided into two groups :

• on the one hand, power shimming group of rods (PSG) that compensates the power defect during load variations, which is the change in core reactivity due to moderator and fuel neutron feedbacks. Four groups of rods compose PSG (fig.

1, left) : grey rods (G1 and G2) and black rods (N1 and N2), which differ by their respective neutron absorption capacity. During a load drop, PSG groups of rods don't plunge into the core simultaneously but with constant overlaps to prevent too much axial power imbalance. Let r 1 (resp. r 2 ) be the overlap value between groups G1 and G2 (resp. G2 and N1) : G2 moves as soon as G1 has reached an insertion step value of r 1 as well as N1 moves when G2 exceeds r 2 inserted steps. The position of the PSG is slaved to the net electrical power of the nuclear power plant ;

• on the other hand, temperature regulation group of rods (TRG) aims at respecting the temperature program set by the operating technical specifications all along the load transient. NPP operator inserts (resp. extracts) the TRG as soon as core water is hotter (resp. colder) than the setpoint with a dead band of 0.8°C. The operator can move TRG rods all along a fixed manoeuvring band whose width is given in steps by the parameter b, 27 steps for a typical P4 1300 MW reactor [START_REF] Grard | Physique, fonctionnement et sûreté des REP : Le réacteur en production[END_REF].

Boron

Soluble boron is a neutron poison diluted in core's primary water. Contrary to the rods its concentration adjustment is much slower than a rod step movement and its effect are homogeneous in the core. In G mode load follow control, it is mainly used to compensate reactivity variations occurring during the evolution of xenon concentration in the core, which is a strong neutron absorber with a delayed dynamic at the time scale of a load drop. Moreover, xenon equilibrium of production (mostly due to iodine decay) and consumption (by neutron capture) is highly tied to the variation of neutron flux in the core. Not only a power drop diminishes neutron flux absolutely, it also redistributes core's axial power profile due to combined effects of rods insertion and an axially evolutionary negative moderator temperature coefficient. In case of reactivity loss, the plant operator injects clear water to dilute boron whereas he increases boron concentration in case of abnormal reactivity rise.

Control of Load Variations

Core Computer Simulator in Load-Following Mode

G mode load control consists in the coordinated use of both removable neutron absorbing rods and soluble boron diluted in the water of the primary circuit. Basically, controlling a nuclear reactor consists in adjusting the position of the rods and the boron concentra-tion to enable load variations while respecting safety margins at all times. Much effort has been put into simulating the operation of a nuclear reactor. Strategies aiming at modeling a human behavior regarding the interpretation of the observables of the core and the decisions which result from it have been developed [START_REF] Muniglia | A multi-physics pwr model for the load following[END_REF]Drouet et al., 2020a).

Similarly, PID controllers1 performs well to optimize the current control mode and minimize the power oscillations mainly related to the axial transients of the xenon [START_REF] Zeng | A functional variable universe fuzzy pid controller for load following operation of pwr with the multiple model[END_REF][START_REF] Rafiei | Load-following control of a nuclear reactor using optimized fopid controller based on the twopoint fractional neutron kinetics model considering reactivity feedback effects[END_REF][START_REF] Mahfudin | Auto-tuning pid controller for nuscale nuclear reactor using point reactor kinetics model simulator[END_REF]. In this study, we choose to opt for a model close to the real case G mode and we search for an optimization focusing on variables that do not intrinsically affect the behavior of the dummy operator for a corrective action using rods or boron. Whatever the model of a fictitious operator, safety rules impose to guarantee the axial power stability of the core during load following mode. The observable parameter ∆I quantifies it ; given the power produced in the core's upper half and lower half denoted P U and P D respectively, and P rel the relative power of the core :

∆I = P U -P D P U + P D P rel , given in % of core's height. (1)
A plant operator must control the reactor during load following so that ∆I remains within a margin of ± 5% of its value at its state before the power drop. Maintaining a good axial power stability may lead the operator to favor control of the neutron flux by the boron rather than by the rods, since they distort the power profile of the core more widely. However, relying more extensively on boron adjustments in the primary circuit means a rise in fluid transfers leading to the generation of more liquid waste one seeks to minimize for technical, economic and environmental concerns.

Optimization Problem Setting

Previous work has introduced solutions of better sets of PSG overlaps and width of the TRG manoeuvring band for a 6-6 hours transient (6 hours at 30% of nominal power followed by 6 hours at full power) to guarantee a better stability against axial displacement of the neutron flux distribution as well as on the quantity of liquid effluents produced during the load transient [START_REF] Muniglia | A fitness landscape view on the tuning of an asynchronous master-worker ea for nuclear reactor design[END_REF]Drouet et al., 2020a). Finding the best set of operating levers whatever the amplitude of load drops is a perspective of these studies.

We introduce in this work the Nuclear Unit Commitment Optimization Problem (NUCOP), which aims to jointly find the best distribution of a fleet load transient to a group of two reactors while optimizing rods parameters. Solving NUCOP would lead to judge the better way of meeting grid's needs with a given fleet of reactors while staying within safety rules and minimizing both objectives as it investigates two main issues :

• an optimized way to involve both reactors to assume the fleet transient instead of letting a single reactor doing a deep load drop (fig. 2) ;

• the joint search of the best division of the fleet transient with the best set of PSG rods overlaps and TRG manoeuvring band's width adapted for both reactors. Search space : NUCOP is a discrete and expensive bi-objective optimization problem with a 5dimensional search space of size ∼ 10 10 , conducted on a fleet composed of two French PWR reactors of 1300 MW each. First type of search space's variables are rods overlaps r 1 and r 2 as well as the width b of the maneuvering band of TRG. These numerical variables are counted into rod steps. Second type of variables are categorical and identify the unit commitment dispatch of the two-reactor fleet. As shown on figure 2, we enable each units to operate any load drop as long as the whole fleet drop reaches 50% of fleet's nominal power. The fleet drop is divided into two time zones which leads to two variables of power distribution d 1 and d 2 . For each d i we fix 40 possible possibilities for the power distribution. Search space is summarized in table 1.

Fitness space : Given the domain X = {(x 0 , . . . , x 4 ), x i ∈ E i } ⊂ N 5 , the objective is to find all Pareto-optimal solutions x * such as :

x * = arg min x∈X f (x)
(2)

E 0 = [255] x 1 r 2 E 1 = [255] x 2 b E 2 = [7, 118] x 3 d 1 E 3 = [39] x 4 d 2 E 4 = [39]
where f : X → R 2 is a black-box function computationally expensive to evaluate which returns a two objective values. Objective functions respectively stands for the Euclidean distances in the 2D space of the generated liquid effluents (V 1 ,V 2 ), and the integral of axial power imbalance (∆I 1 , ∆I 2 ) during each transient and for each reactor. The black-box fitness function is :

f : X → R 2 x → z = f 1 (x) f 2 (x) =   V 1 2 +V 2 2 ∆I 1 2 +∆I 2 2   .
(3)

4 Optimization using Adaptive MOEA/D with Walsh Surrogate

The NUCOP problem is a bi-objective black-box optimization problem. This kind of problem can be solved with a variety of metaheuristic methods whose goal is to find the set of non-dominated solutions of the problem. For a bi-objective minimization problem, a solution x is dominated by

x ′ iff f 1 (x ′ ) ≤ f 1 (x) and f 2 (x ′ ) ≤ f 2 (x)
, and if there is an i ∈ {1, 2} such that f i (x ′ ) < f i (x). A solution x ⋆ is non-dominated if there is no other solution x such that x ⋆ is dominated by x. The set of non-dominated solutions is called the Pareto set, the image of this set by f constitutes the Pareto front of the objective space. The metaheuristic used for solving NUCOP is a variant of the MOEA/D algorithm [START_REF] Zhang | Moea/d: A multiobjective evolutionary algorithm based on decomposition[END_REF] which is an evolutionary algorithm based on the decomposition of the objective space into several directions.

MOEA/D for solving NUCOP

Let Λ = {λ 1 , ..., λ Λ } be the set of directions splitting the objective space. MOEA/D consists in dividing the bi-objective problem into Λ single objective subproblems communicating with each other to enable a cooperatively solving of the initial problem. Each subproblem aims at minimizing f along its respective assigned direction. The single objective function associated to each direction λ q is the Chebyshev function defined as

g(z | λ q ) = max i∈{1,2} λ q,i |z ref i -z i | (4)
where z ref is the reference point set to the minimal values of objectives. We save in a buffer file each best solutions found for each subproblem (namely for each direction λ q ). Each iteration of the algorithm (see Algorithm 1) gives the fitness values of a candidate solution x in a direction λ q ; we compute g(x | λ q ) and compare the result to solutions saved in the buffer file whose directions are in the set T λ q (T ), which refers to the T neighbouring directions of λ q . We replace any solution of the buffer file whose corresponding direction is in T λ q (T ) by x iff x dominates it.

Algorithm 1: MOEA/D : master process In addition, our MOEA/D algorithm benefits from a massive parallel architecture and follows a masterworkers procedure (algorithms 1 and 2). The computation time to evaluate a candidate solution of NU-COP is expensive (on average 45 minutes). Moreover, this computation time has a large variance according to the rod parameters and fleet's transient distribution.

1 for p ← 1 to N proc do 2 x p ← Init.
Because of limited computing resources, the masterworkers algorithm is asynchronous to avoid wasting time between the evaluation of the first and the last worker associated to each direction.

Algorithm 2: MOEA/D: worker process 1 Receive x p from master process ; 2 Compute z p = f (x p ) ; 3 Send (x p , z p ) to master process ;

Walsh Surrogate Model

One way to face expensive fitness evaluation is to use a surrogate model in order to extract more information from the costly black-box fitness function, and to better guide the search. Popular methods of assisted optimization with surrogate models try to substitute the original time-consuming black-box function with a cheap approximation of it : the best solution given by the surrogate function is then computed as the next candidate through the real problem function. In the context of discrete optimization, Walsh functions have already been used as surrogate model for solving pseudo-boolean optimization problems [START_REF] Leprêtre | Walsh functions as surrogate model for pseudoboolean optimization problems[END_REF][START_REF] Verel | A surrogate model based on walsh decomposition for pseudo-boolean functions[END_REF].

Walsh functions

Notations : we reuse range notations described in table 1. Also, ♯E refers to the cardinality of a set E.

The set of Walsh functions (W k ) k∈[2 d -1] is a normal and orthogonal basis of the space of pseudo-Boolean functions F d = { f : {0, 1} d → R} with d being an integer [START_REF] Walsh | A closed set of normal orthogonal functions[END_REF][START_REF] Bethke | Genetic algorithms as function optimizers[END_REF], [START_REF] Leprêtre | Walsh functions as surrogate model for pseudoboolean optimization problems[END_REF][START_REF] Verel | A surrogate model based on walsh decomposition for pseudo-boolean functions[END_REF]. Expansions based on Walsh functions are sometime called Fourier expansions as they have similar [START_REF] O'donnell | Analysis of boolean functions[END_REF]. By virtue of basis properties, any f ∈ F d can be written as a unique linear combination of Walsh functions:

∀x ∈ {0, 1} d , f (x) = 2 d -1 ∑ k=0 ω k W k (x) (5)
where ω k ∈ R is the coefficient of the k-th Walsh function W k defined by:

W k (x) = (-1) ∑ d-1 i=0 k i x i (6)
where k i ∈ {0, 1} is the i th binary digit of the integer k. As classical Fourier series for an analog signal, and in order to reduce the number of terms in the series, it is possible to truncate the Walsh expansion of f to get a polynomial approximation function f . We call order of the k-th Walsh function the number of 1 in the binary representation of k. The Walsh function of order 0 is the constant term, Walsh functions of order 1 are linear functions as they only depend on 1 binary variable, Walsh functions of order 2 are quadratic, and so on. This leads to the following rearrangement of the series:

f (x) = 2 d -1 ∑ k=0 ω k W k (x) = d ∑ n=0 ∑ k with ord(W k )=n ω k W k (x). (7)
It should be noted that Walsh's surrogate model can be viewed as a linear regression in a specific Hilbertian space.

One-Hot Encoded Boolean Representation of the Search Space

We consider a one-hot encoded Boolean representation of any variable of x = (x 0 , x 1 , x 2 , x 3 , x 4 ) using an uniform partitioning of each range E i corresponding to the domain of the variable x i (see tab. 1). We define for each variable x i a binary digit sequence of size ♯E i -1 (see eq. 8) such that at most one digit is equal to 1. More formally, let be (E i j ) an uniform partition of E i into ♯E i sets of equal cardinality : E i = ♯E i j=0 E i j , and E i j 1 ∩E i j 2 = / 0 for any j 1 ̸ = j 2 . The proposed onehot encoded strategy follows two rules (eq. 9). Firstly the unique digit sequence with no 1 corresponds to E i 's first element2 , secondly other elements are linked to sequences that have exactly one 1 ; its position is given by the value x i : if x i ∈ E i j for j ∈ [1, ♯E i ], then the digit 1 is positioned at index j -1 in the boolean sequence. Other digits equal to 0.

h i : E i → {0, 1} ♯E i x i → x i 0 , ..., x i(♯E i -1) (8) ∀ j ∈ [♯E i -1], x i j = 1 if x i ∈ E i( j+1) 0 otherwise. ( 9
)
The one-hot encoded Boolean representation of the whole variable x = (x 1 , ..., x 5 ) results from the concatenation of one-hot encoded sequences corresponding to each variable x i : h(x) = (h 0 (x 0 ), h 1 (x 1 ), . . . , h 4 (x 4 )).

Example : let x = (185, 175, 27, 4, 36) a random element of X . If ♯E i = 4 for all i then for instance the uniform partitioning of r 1 pa- 192,255], which means that h(x 0 ) = (0010). Following this schemata with other subvariables the complete Boolean sequence of x would be (0010 0010 1000 1000 0001).

rameter domain is E 0 = [255] = [63] ∪ [64, 127] ∪ [128, 191] ∪ [
Such a one-hot encoded representation is common with previous work proposed recently in [START_REF] Dadkhahi | Fourier representations for black-box optimization over categorical variables[END_REF], used for Boolean Fourier representation of any categorical variable. It is very practical as we can remove many terms in the expression of f (x) (eq. 7) ; indeed, for any sub-variable x i , Walsh terms or order n ≥ 2 involving more than two elements of the same family (x i j ) are inconsistent. Considering that NUCOP is a d = 5 dimension problem (tab. 1) :

f (x)=ω 0 + d ∑ n=1 ∑ I∈ [1, d] n ∑ × ♯I c=1 [♯EI c -1] ω I,J W I,J (x) (10)
where [1, d] n

is
W I ,J (x) = ∏ (i, j)=(I c ,J c ) with c∈ 1,#I (-1) x i j . ( 11
)
As mentioned in equation ( 7) we can sort the terms of the series (10) according to their contribution in linear or interaction effects within the linear regression equation that defines Walsh surrogate model. Firstorder terms of ( 10) are all those whose Walsh function W I ,J (x) involves a single x i j ; second-order terms are terms containing two terms in the product (11) ; and so forth. We can write the truncated expansion f (x) of f (x) at order n = 2 by the Walsh decomposition given in (10) :

f (x) = ω0 + ∑ i∈[1,5] ∑ j∈[♯E i -1] ωi j (-1) x i j + ∑ (i,i ′ )∈ 1, 5 n ∑ ( j, j ′ )∈[♯E i -1]×[♯E i ′-1] ωi j,i ′ j ′ (-1) x i j +x i ′ j ′ (12)
Example : consider a Walsh decomposition model based on x 4 and x 5 sub-variables of the NUCOP search space. Let's partition E 4 and E 5 with cardinality ♯E 4 = 2 and ♯E 5 = 3 respectively. Approxima-tion (12) of f (x) gives :

f (x) = ω0 + ω4,0 (-1) x 4,0 + ω4,1 (-1) x 4,1 + ω5,0 (-1) x 5,0 + ω5,1 (-1) x 5,1 + ω5,2 (-1) x 5,2 + ω4,0,5,0 (-1) x 4,0 +x 5,0 + ω4,0,5,1 (-1) x 4,0 +x 5,1 + ω4,0,5,2 (-1) x 4,0 +x 5,2 + ω4,1,5,0 (-1) x 4,1 +x 5,0 + ω4,1,5,1 (-1) x 4,1 +x 5,1 + ω4,1,5,2 (-1) x 4,1 +x 5,2 . (

Following semantics proposed in [START_REF] Dadkhahi | Fourier representations for black-box optimization over categorical variables[END_REF], we refer to the tuple (x i j ) j∈[♯E i -1] as the descendant of the parent variable x i of the solution x ∈ X . The completeness and uniqueness of the decomposition ( 10) is justified by the fact that each W I ,J (x) contains exactly one x i j descending from the parent variable x i . The proof is given in Appendix of [START_REF] Dadkhahi | Fourier representations for black-box optimization over categorical variables[END_REF].

Walsh Adaptive Method

As MOEA/D is an evolutionary algorithm, the mutation step involves the use of a selection and a mutation operator. Given a evaluated solution x and a mutation rate vector P = (P 1 , ..., P 5 ) T , the first one performs a Bernoulli test with an expectancy value of P i on each variable x i of x, so that sub-variable x i is mutated into xi iff the test succeeds. Usually, P is fixed such as each sub-variable has the same probability to mutate, which means that each P i equals to 1/d = 0.2. The second one is based on a vector R = (R 1 , ..., R 5 ) T such that all mutated variables xi are chosen thanks to an uniform distribution in the real interval [x i -R i ⌊ub ilb i ⌋, x i + R i ⌊ub ilb i ⌋], with ub i and lb i respectively being the upper bound and lower bound values of x i .

Our adaptive method is based on an evolutionary selection operator where vector P is regularly updated during a run, following the evolution of the Walsh surrogate model coefficients. In this study, P updates occur once all worker processors deployed for a run of NUCOP resolution have returned a solution to the master processor. The adaptive procedure starts by computing the current Pareto optimal solutions that will be used as the sample base for training the surrogate model. Then we use python library Scikit-Learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] to build a surrogate for each one of the objectives under study, which consists in a linear regression fit using Lasso interpolator. As mentioned in [START_REF] Leprêtre | Walsh functions as surrogate model for pseudoboolean optimization problems[END_REF] the use of Lasso is motivated by the high number of Walsh coefficients in the decomposition, as it relies on the choice of the partitioning mesh of each E i set. For instance, if eight sets make a partition of the interval E 0 = [255],

there will be a family ( ω1, j ) j∈ [7] of eight coefficients associated with parameter r 1 in the Walsh surrogate model. After fitting the surrogate we get access to all Walsh coefficients that can be separated according to the sub-variable they correspond with : taking the L1-norm on coefficients linked to the same x i , we perform a normalization by the L1-norm of all Walsh surrogate coefficients (except for ω 0 which a constant term). This reasoning leads to the calculation of two new selection vectors Pf 1 and Pf 2 respectively for objective f 1 and f 2 , that involves a learning rate coefficients denoted l. For instance, in the case of objective f 1 :

P f 1 = (1 -l)P f 1 ,i + l Pf 1 ,i , ∀i ∈ [4] . (14)
The learning rate parameter influence the importance given to the modification of P f 1 during iterations. In this study, we have taken l = 0.2 for all runs because of limited computational resources, but with more computational budget it would possible to carry out a parametric study on this parameter in order to assess its effects on the study proposed here. We use Pf 1 and Pf 2 to update regularly the mutation rate vector of MOEA/D such as it is given by :

P = (1 -l)P i + l Pf 1 ,i + Pf 2 ,i 2 , ∀i ∈ [4] . (15)
The underlying idea is to progressively favor the selection of variables whose mutation is expected to have more pronounced effects on the results given by the black-box function. In short, we take advantage of the interpretation of the coefficients which impact more the objective values in a linear regression to adapt the selection operator of MOEA/D (see alg.

3).

Algorithm 3: Adaptive mutation operator

EXPERIMENTAL ANALYSIS

Experimental setup

The experimental setup is as follows. We use 600 computational processors, distributed over 200 directions decomposing the objective space (3 per direction). The neighborhood of a direction is equal to the total number of directions in the problem to enable a total communication between directions following recommendations of previous works working on the optimization of rod parameters [START_REF] Muniglia | A fitness landscape view on the tuning of an asynchronous master-worker ea for nuclear reactor design[END_REF][START_REF] Drouet | Surrogateassisted asynchronous multiobjective algorithm for nuclear power plant operations[END_REF]. We use Chebyshev's scalarization function to evaluate the dominance of one solution over another in the objective space. For each raw fit from the controlling simulator, we center-reduce the fit according to a normal distribution N (µ, σ) whose parameters are derived from a totally random sampling of the search space upstream of the NUCOP optimization study. The proposed adaptive mutation operator method is denoted AMOW for adaptive mutation operator based on Walsh decomposition, and it compared the baseline method FMO for fixed mutation operator. The AMOW method updates the mutation operator every M = 600 evaluated solutions. We train a Walsh surrogate model by evaluating current non-dominated solutions that are taken to build a training dataset with the corresponding fitness values as training targets. Each sample in the dataset is then expressed as a feature of five features (r 1 , r 2 , b, d 1 and d 2 ) and two targets ( f 1 and f 2 ). The construction of the Walsh predictors made of all W I ,J leads to a nonorthogonal and non-normal basis : as shown on figure 3 the non-transformation of the base vectors through the Gram Schmidt algorithm still slightly decreases the mean absolute error observed on both objectives ; thus, results presented in the rest of this study were obtained with the non-orthonormal basis of Walsh functions3 . It should be noted that Walsh surrogate has similar performance as decision tree method ; in any case such a method is rejected as it does not give access to the individual contribution of each research space's sub-variable.

We perform 10 runs that computes 9000 solutions each. The performance assessment protocol is such that we compare two by two runs that begin with the same initial population given by an latin hypercube sampling of the search space (see table 1 for dimensions). In other words, there are 10 different runs for each method FMO and method AMOW, and we use only ten different population samples to be-gin the each run. Performance indicators are given as the mean value observed on all runs. Concerning the NUCOP problem, we study a fleet of two reactors : one reactor has got fresh fuel whereas the other one is at the middle of its fuel campaign. Calculations are performed at the TGCC on Intel Skylake AVX-512 2.7GHz CPU units (GENCI, 2023).

Performance Analysis

Figure 4 shows the evolution of the coordinates of the mutation rate vector P during the run. The optimization process progressively reduces the probability that b and r 2 sub-variables will be mutated for the next 600 iterations, thus favoring mutation of the other subvariables. In other words, we take advantage of a quantification of the linear effect induced by variation of b on objectives thanks to Walsh's model that privilege mutation of the other parameters with the hope of exploring the fitness landscape more effectively. In the field of nuclear reactor physics for power maneuvering, this result shows that there is less expected improvement in searching for the best width of the maneuvering strip of TRG rods than in changing other variables like rod overlaps and fleet transient's distribution. The physical interpretation of these interactions belongs to the field of reactor physics that are beyond the scope of this current study.

One notable result is that decreasing the chances of mutating variable b does not seem to have much impact on the diversity of solutions in the compromise zone between the two objectives (figure 6). We believe the use of AMOW makes it possible to progressively remove a variable from the NUCOP problem if it turns out that its mutations are not likely to cause a large variation of the black-box function, without leading to a great loss of potential non-dominated solutions. In the formalism of linear regression, this amounts to discarding the terms with the weakest main effects. Figure 6 also shows that best solutions minimizing absolutely f 1 or f 2 were done with the AMOW method. However these results must be qualified to the extent that most of the metrics usually used for the comparison of optimization method seems to favour the FMO non-adaptive method (hypervolume metric is presented on figure 5). However, a Welsch ttest made on the hypervolume metric gives a p-value of 0.2606 meaning that we cannot statistically state that either method leads to better performance on the hypervolume metric. The use of AMOW seems more relevant to get a better comprehension of the physics of the fitness function by reducing the search space than for hoping to achieve a better approximation of the optimal Pareto front. Points stand for mean values over all runs whereas shaded areas stand for the 1σ confidence interval (note that y-axis scale is logarithmic so mean values are not centered in error bars).

CONCLUSION AND PERSPECTIVES

Modeling a nuclear reactor is complex and simulating load following operations is costly in terms of computing time. The joint optimization of the control rod overlap settings, TRG maneuvering bandwidth and the distribution of a fleet transient over a group of two reactors is an innovative problem whose resolution leads to new considerations on the possible gains in nuclear power maneuverability. In this context, we have implemented within a MOEA/D algorithm an adaptive mutation rate based on a decomposition of the fitness function in the basis of Walsh functions. This method favours the mutation of variables most likely to lead to greater goal variation at the expense of those contributing least to the problem, as if a fitness landscape analysis is made in-line instead of being done before the optimization process. Naturally, the adaptive method using the surrogate model is not limited to load following operation in nuclear energy context. Future works will test the proposed adaptive method to others discrete expensive problems. Future research works may involve developing an adaptive mutation operator up to the second order of the Walsh series decomposition, to take advantage of knowledge of the interactions between variables in the problem, and moreover mutation operators which mutate several variables could also be designed in order to increase the convergence rate of the optimization algorithms for larger dimension problems.

Figure 1 :

 1 Figure 1: Typical load map for 1300 MW French PWR (left) and zoom on a single assembly provided with black rods (right).

Figure 2 :

 2 Figure 2: Example of a two-reactor distribution for the NU-COP problem. Fleet transient on the left is divided into two transients on the right. In this case, both reactors have the same nominal power. Amplitudes are not to scale.

  the set of subsets of cardinal n drawn from the set [1, d] and × ♯I c=1 [♯E I c -1] is the cartesian product of all integer sets [♯E I c -1] with I c taking values in [1, d] according to the sum term exploring each n-subset previously defined. Finally, Walsh functions W I ,J are given by

Figure 3 :

 3 Figure 3: Mean absolute error on f 1 (top) and f 2 (bottom) obtained from different surrogate methods during a run. The training dataset is updated each 600 tested solutions.Points stand for mean values over all runs whereas shaded areas stand for the 1σ confidence interval (note that y-axis scale is logarithmic so mean values are not centered in error bars).

Figure 4 :Figure 5 :

 45 Figure4: Evolution of mutation rate vector's coordinates during the run, using method AMOW. Interpretations of points and shaded areas are the same as in figure3.

Figure 6 :

 6 Figure6: Pareto fronts for all runs (method FMO in blue, method AMOW in red).

Table 1 :

 1 Search space sub-variables and their range.Brackets [a, b] refers to the integer range {a, a + 1, ..., b} while shorter notation [b] is for the {0, 1, ..., b} range.

	Sub-variable Parameter	Range
	x 0	r 1
	.

PID controllers computes the corrective action by estimating the (p)roportionnal, (i)ntegral and (d)erivative signal error over time.

All sets E i are sorted in ascending order so that considering the first element of E i makes sense.

Errata (underlined words) : study shows results using Walsh original basis, not the orthonormal one. Still, experiments using both methods have led to similar results.

ACKNOWLEDGEMENTS

Authors would like to thank the GENCI project for access to the CCRT's computing resources, which were used during allocation campaign A0131013043 (GENCI, 2023).