Christophe Alias

Improving Communication Patterns in Polyhedral Process Networks

Keywords: High-level synthesis, polyhedral compilation, FIFO, FPGA Synthèse de circuit haut-niveau, compilation polyédrique, FIFO, FPGA

Embedded systems performances are bounded by power consumption. The trend is to offload greedy computations on hardware accelerators as GPU, Xeon Phi or FPGA. FPGA chips combine both flexibility of programmable chips and energy-efficiency of specialized hardware and appear as a natural solution. Hardware design is long, fastidious and bug prone. Hardware compilers from high-level languages (High-level synthesis, HLS) are required to exploit all the capabilities of FPGA while satisfying tight time-to-market constraints. Compiler optimizations for parallelism and data locality restructure deeply the execution order of the processes, hence the read/write patterns in communication channels. This breaks most FIFO channels, which have to be implemented with addressable buffers. Expensive hardware is required to enforce synchronizations, which often results in dramatic performance loss. In this paper, we present an algorithm to partition the communications so that most FIFO channels can be recovered after a loop tiling, a key optimization for parallelism and data locality. Experimental results show a drastic improvement of FIFO detection for regular kernels at the cost of (few) additional storage. As a bonus, the storage can even be reduced in some cases.

Introduction

Since the end of Dennard scaling, the performance of embedded systems is bounded by power consumption. The trend is to trade genericity (processors) for energy efficiency (hardware accelerators) by offloading critical tasks to specialized hardware. FPGA chips combine both flexibility of programmable chips and energy-efficiency of specialized hardware and appear as a natural solution. High-level synthesis (HLS) techniques are required to exploit all the capabilities of FPGA, while satisfying tight time-to-market constraints. Parallelization techniques from high-performance compilers are progressively migrating to HLS, particularly the models and algorithms from the polyhedral model -a powerful framework to design compiler optimizations. Additional constraints must be fulfilled before plugging a compiler optimization into an HLS tool. Unlike software, the hardware size is bounded by the available silicum surface. The bigger is a parallel unit, the less it can be duplicated, thereby limiting the overall performance. In particular, tough program optimizations are likely to spoil the performances if the circuit is not post-optimized carefully [START_REF] Alias | Optimizing Affine Control with Semantic Factorizations[END_REF]. An important consequence is that the the roofline model is not longer valid in HLS [START_REF] Bruno Da Silva | Performance modeling for fpgas: extending the roofline model with high-level synthesis tools[END_REF]. Indeed, peak performance is no longer a constant: it decreases with the operational intensity. The bigger is the operational intensity, the bigger is buffer size and the less is the space remaining for the computation itself. Consequently, it is important to produce at source-level a precise model of the circuit which allows to predict accurately the resource consumption. Process networks are a natural and convenient intermediate representation for HLS [START_REF] Alias | Data-aware Process Networks[END_REF][START_REF] Turjan | Compiling nested loop programs to process networks[END_REF][START_REF] Verdoolaege | Polyhedral Process Networks[END_REF]. A sequential program is translated to a process network by partitioning computations into processes and flow dependences into channels. Then, the the process and buffers are factorized and mapped to hardware.

In this paper, we focus on the translation of buffers to hardware. We propose an algorithm to restructure the buffers so they can be mapped to inexpensive FIFOs. Most often, a direct translation of a regular kernel -without optimization -produces to a process network with FIFO buffers [START_REF] Turjan | Classifying interprocess communication in process network representation of nested-loop programs[END_REF]. Unfortunately, data transfers optimization [START_REF] Alias | Optimizing remote accesses for offloaded kernels: Application to high-level synthesis for FPGA[END_REF] and generally loop tiling reorganizes deeply the computations, hence the read/write order in channels (communication patterns). Consequently, most channels may no longer be implemented by a FIFO. Additional circuitry is required to enforce synchronizations [START_REF] Alias | Data-aware Process Networks[END_REF][START_REF] Zissulescu | Solving out of order communication using cam memory; an implementation[END_REF][START_REF] Turjan | Realizations of the extended linearization model. Domain-specific processors: systems, architectures, modeling, and simulation[END_REF][START_REF] Van Haastregt | Enabling automatic pipeline utilization improvement in polyhedral process network implementations[END_REF] which result in larger circuits and causes performance penalties. In this paper, we make the following contributions:

We propose an algorithm to reorganize the communications between processes so that more channels can be implemented as FIFO after a loop tiling. As far as we know, this is the first algorithm to recover FIFO communication patterns after a compiler optimization.

Experimental results show that we can recover most of the FIFO disabled by communication optimization, and more generally any loop tiling, at almost no extra storage cost.

The remainder of this paper is structured as follows. Section 2 introduces polyhedral process network and discusses how communication patterns are impacted by loop tiling, Section 3 describes our algorithm to reorganize channels, Section 4 presents experimental results, Finally, Section 5 concludes this paper and draws future research directions. 1

Preliminaries

This section defines the notions used in the remainder of this paper. Section 2.1 and 2.2 introduces the basics of compiler optimization in the polyhedral model and defines loop tiling. Section 2.3

Alias

defines polyhedral process networks (PPN), shows how loop tiling disables FIFO communication patterns and outlines a solution.

Polyhedral Model at a Glance

Translating a program to a process network requires to split the computation into processes and flow dependences into channels. The polyhedral model focuses on kernels whose computation and flow dependences can be predicted, represented and explored at compile-time. Specifically, the control must be predictable: (for loops, if with conditions on loop counters), arrays only; and loop bounds, conditions and array accesses must be affine functions of surrounding loop counters and structure parameters (typically the array size). This way, the computation may be represented with Presburger sets (typically approximated with convex polyhedra, hence the name). This makes possible to reason geometrically about the computation and to produce precise compiler analysis thanks to integer linear programming: flow dependence analysis [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF], scheduling [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] or code generation [START_REF] Bastoul | Efficient code generation for automatic parallelization and optimization[END_REF][START_REF] Quilleré | Generation of efficient nested loops from polyhedra[END_REF] to quote a few. Most compute-intensive kernels from linear algebra and image processing fit in this category. In some case, kernels with dynamic control can even fit in the polyhedral model after a proper abstraction [START_REF] Alias | Multi-dimensional rankings, program termination, and complexity bounds of flowchart programs[END_REF]. Figure 1.(a) depicts a polyhedral kernel and (b) depicts the geometric representation of the computation for each assignment (• for assignment load, • for assignment compute and • for assignment store). The vector i = (i 1 , . . . , i n ) of loop counters surrounding an assignment S is called an iteration of S.

The execution of S at iteration i is denoted by S, i . The set D S of iterations of S is called iteration domain of S. The original execution of the iterations of S follows the lexicographic order over D S . For instance, on the statement C: (t, i) (t , i ) iff t < t or (t = t and i < i ). The lexicographic order over Z d is naturally partitioned by depth: = 1 . . .

d where (u 1 . . . u d ) k (v 1 , . . . , v d ) iff ∧ k-1 i=1 u i = v i ∧ u k < v k .
Dataflow Analysis On Figure 1.(b), red arrows depict several flow dependences (read after write) between executions instances. We are interested in flow dependences relating the production of a value to its consumption, not only a write followed by a read on the same location. These flow dependences are called direct dependences. Direct dependences represent the communication of values between two computations and allow to rule communications and synchronizations in the final process network. They are crucial to build the process network. Direct dependences can be computed exactly in the polyhedral model [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF]. The result is a relation → relating each producer P, i to one or more consumers C, j . Technically, → is a Presburger relation between vectors (P, i) and vectors (C, j) where assignments P and C are encoded as integers. For example, dependence 5 is summed up with the Presburger relation:

{(t -1, i) → (t, i), 0 < t ≤ T ∧ 0 ≤ i ≤ N }.
Presburger relations are computable and efficient libraries allow to manipulate them [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF][START_REF] Kelly | The omega calculator and library, version 1.1. 0[END_REF]. In the remainder, direct dependence will be referred as flow dependence or dependence to simplify the presentation.

Scheduling and Loop Tiling

Compiler optimizations change the execution order to fulfill multiple goals such as increasing the parallelism degree or minimizing the communications. The new execution order is specified by a schedule. A schedule θ S maps each execution S, i to a timestamp θ S ( i) = (t 1 , . . . , t d ) ∈ Z d , the timestamps being ordered by the lexicographic order . In a way, a schedule dispatches each execution instance S, i into a new loop nest, θ S ( i) = (t 1 , . . . , t d ) being the new iteration vector of S, i . A schedule θ induces a new execution order ≺ θ such that S, i ≺ θ T, j iff θ S ( i) θ T ( j). Also, S, i θ T, j means that either S, i ≺ θ T, j or θ S ( i) = θ T ( j). When a schedule is injective, it is said to be sequential: no execution is scheduled at the same time, hence everything is executed in sequence. In the polyhedral model, schedules are affine functions. They can be derived automatically from flow dependences [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. On Figure 1, the original execution order is specified by the schedule θ load (i) = (0, i), θ compute (t, i) = (1, t, i) and θ store (i) = (2, i). The lexicographic order ensures the execution of all the load instances (0), then all the compute instances ( 1) and finally all the store instances (2). Then, for each statement, the loops are executed in the specified order. Loop tiling is a transformation which partitions the computation in tiles, each tile being executed atomically. Communication minimization [START_REF] Alias | Optimizing remote accesses for offloaded kernels: Application to high-level synthesis for FPGA[END_REF] typically relies on loop tiling to tune the ratio computation/communication of the program beyond the ratio peak performance/communication bandwidth of the target architecture. Figure 2.(a) depicts the iteration domain of compute and the new execution order after tiling loops t and i. For presentation reasons, we depict a domain bigger than in Figure 1.(b) (with bigger N and M ) and we depict only a part of the domain. In the polyhedral model, a loop tiling is specified by hyperplans with linearly independent normal vectors τ 1 , . . . , τ d where d is the number of nested loops (here τ 1 = (0, 1) for the vertical hyperplanes and τ 2 = (1, 1) for the diagonal hyperplanes). Roughly, hyper-plans along each normal vector τ i are placed at regular intervals b i (here b 1 = b 2 = 2) to cut the iteration domain in tiles. Then, each tile is identified by an iteration vector (φ 1 , . . . , φ d ), φ k being the slice number of an iteration i along normal vector τ k :

Inria for i := 0 to N + 1 • load(a[0, i]); for t := 1 to T for i := 1 to N • a[t, i] := a[t -1, i -1] + a[t -1, i]+ a[t -1, i + 1]; for i := 1 to N • store(a[T, i]);
φ k = τ k • i ÷ b k . The result is a new iteration domain (here D = {(φ 1 , φ 2 , t, i), 2φ 1 ≤ t < 2(φ 1 + 1) ∧ 2φ 2 ≤ t + i < 2(φ 2 + 1)}).
In turn, this domain can be processed with polyhedral analysis: the polyhedral model is closed under loop tiling. In particular, the tiled domain can be scheduled. For instance, θS (φ 1 , φ 2 , t, i) = (φ 1 , φ 2 , t, i) specifies the execution order depicted in Figure 2.(a)): tile with point (4,4) is executed, then tile with point (4,8), then tile with point [START_REF] Alias | Data-aware Process Networks[END_REF][START_REF] Quilleré | Generation of efficient nested loops from polyhedra[END_REF], and so on. For each tile, the iterations are executed for each t, then for each i.

Polyhedral Process Networks

Given the iteration domains and the flow dependence relation, →, we derive a polyhedral process network by partitioning iterations domains into process and flow dependence into channels. More formally, a polyhedral process network is a couple (P, C) such that:

Each process P ∈ P is specified by an iteration domain D P and a sequential schedule θ P inducing an execution order ≺ P over D P . Each iteration i ∈ D P realizes the execution instance µ P ( i) in the program. The processes partition the execution instances in the program: {µ P (D P ) for each process P } is a partition of the program computation.

Each channel c ∈ C is specified by a producer process P c ∈ P, a consumer process C c ∈ P and a dataflow relation → c relating each production of a value by P c to its consumption by C c : if i → c j, then execution i of P c produces a value read by execution j of C c . → c is a subset of the flow dependences from P c to C c and the collection of → c for each channel c between two given processes P and C, {→ c , (P, C) = (P c , C c )}, is a partition of flow dependences from P to C.

The goal of this paper is to find out a partition of flow dependences for each producer/consumer couple (P, C), such that most channels from P to C can be realized by a FIFO.

Figure 1.(c) depicts the PPN obtained with the canonical partition of computation: each execution S, i is mapped to process P S and executed at process iteration i: µ P S ( i) = S, i . For presentation reason the compute process is depicted as C. Dependence depicted as i on the dependence graph in (b) are solved by channel i. To read the input values in parallel, we use a different channel per couple producer/read reference, hence this partitioning. We assume that, locally, each process executes instructions in the same order than in the original program: θ load (i) = i, θ compute (t, i) = (t, i) and θ store (i) = i. Remark that the leading constant (0 for load, 1 for compute, 2 for store) has disappeared: the timestamps only define an order local to their process: ≺ load , ≺ compute and ≺ store . The global execution order is driven by the dataflow semantics: an operation is executed whenever its operands are available. The next step is to detect communication patterns to figure out how to implement channels.

Communication Patterns A channel c ∈ C might be implemented by a FIFO iff the consumer C c read the values from c in the same order than the producer P c write them to c (in-order) and each value is read exactly once (unicity) [START_REF] Turjan | Compiling nested loop programs to process networks[END_REF][START_REF] Turjan | Classifying interprocess communication in process network representation of nested-loop programs[END_REF]. The in-order constraint can be written: in-order(→ c , ≺ P , ≺ C ) := ∀x → c x , ∀y → c y : x ≺ C y ⇒ x P y

The unicity constraints can be written:

unicity(→ c ) := ∀x → c x , ∀y → c y : x = y ⇒ x = y
Notice that unicity depends only on the dataflow relation → c , it is independent from the execution order of the producer process ≺ P and the consumer process ≺ C . Furthermore, ¬in-order(→ c , ≺ P , ≺ C ) and ¬unicity(→ c ) amount to check the emptiness of a convex polyhedron, which can be done by most LP solvers. Finally, a channel may be implemented by a FIFO iff it verifies both in-order and unicity constraints:

fifo(→ c , ≺ P , ≺ C ) := in-order(→ c , ≺ P , ≺ C ) ∧ unicity(→ c )
When the consumer reads the data in the same order than they are produced but a datum may be read several times: in-order(→ c , ≺ P , ≺ C ) ∧ ¬unicity(→ c ), the communication pattern is said to be in-order with multiplicity: the channel may be implemented with a FIFO and a register keeping the last read value for multiple reads. However, additional circuitry is required to trigger the write of a new datum in the register [START_REF] Turjan | Compiling nested loop programs to process networks[END_REF]: this implementation is more expensive than a single FIFO. Finally, when we have neither in-order nor unicity: ¬in-order(→ c , ≺ P , ≺ C )∧¬unicity(→ c ), the communication pattern is said to be out-of-order without multiplicity: significant hardware Inria resources are required to enforce flow-and anti-dependences between producer and consumer and additional latencies may limit the overall throughput of the circuit [START_REF] Alias | Data-aware Process Networks[END_REF][START_REF] Zissulescu | Solving out of order communication using cam memory; an implementation[END_REF][START_REF] Turjan | Realizations of the extended linearization model. Domain-specific processors: systems, architectures, modeling, and simulation[END_REF][START_REF] Van Haastregt | Enabling automatic pipeline utilization improvement in polyhedral process network implementations[END_REF]. Consider Figure 1.(c), channel 5, implementing dependence 5 (depicted on (b)) from •, t-1, i (write a[t, i]) to •, t, i (read a[t -1, i]). With the schedule defined above, the data are produced ( •, t -1, i ) and read ( •, t -1, i ) in the same order, and only once: the channel may be implemented as a FIFO. Now, assume that process compute follows the tiled execution order depicted in Figure 2.(a). The execution order now executes tile with point (4,4), then tile with point (4,8), then tile with point [START_REF] Alias | Data-aware Process Networks[END_REF][START_REF] Quilleré | Generation of efficient nested loops from polyhedra[END_REF], and so on. In each tile, the iterations are executed for each t, then for each i. Consider iterations depicted in red as 1, 2, 3, 4 in Figure 2.(b). With the new execution order, we execute successively 1,2,4,3, whereas an in-order pattern would have required 1,2,3,4. Consequently, channel 5 is no longer a FIFO. The same hold for channel 4 and 6. Now, the point is to partition dependence 5 and others so FIFO communication pattern hold. Dependence 5 is partitioned in 3 parts: red dependences crossing tiling hyperplane φ 1 (direction t), blue dependences crossing tiling hyperplane t + i (direction t + i) and green dependences inside a tile. Since the execution order in a tile is the same than the original execution order (actually a subset of the original execution order), green dependences will clearly verify the FIFO communication pattern. As concerns blue and red dependences, source and target are executed in the same order because the execution order is the same for each tile and dependence 5 happens to be short enough. In practice, this partitioning is effective to reveal FIFO channels. In the next section, we propose an algorithm to find such a partitioning.

Our Algorithm

Figure 3 depicts our algorithm for partitioning channels given a polyhedral process network (P, C) (line 5). For each channel c from a producer P = P c to a consumer C = C c , the channel is partitioned by depth along the lines described in the previous section (line 7). D P and D C are assumed to be tiled with the same number of hyperplanes. P and C are assumed to share a schedule with the shape: θ(φ 1 , . . . , φ n , i) = (φ 1 , . . . , φ n , i). This case arise frequently with tiling schemes for I/O optimization [START_REF] Alias | Data-aware Process Networks[END_REF]. If not, the next channel → c is considered (line 6). The split is realized by procedure split (lines 1-4). A new partition is build starting from the empty set. For each depth (hyperplane) of the tiling, the dependences crossing that hyperplane are filtered . θ P (x) ≈ n θ C (y) means that the n first dimensions of θ P (x) and θ C (y) (tiling coordinates (φ 1 , . . . , φ n )) are the same: x and y belong to the same tile. Consider the PPN depicted in Figure 1.(c). with the tiling and schedule discussed above : process compute is tiled as depicted in Figure 2.(c) with the schedule θ compute (φ 1 , φ 2 , t, i) = (φ 1 , φ 2 , t, i). Since processes load and store are not tiled, the only channels processed by our algorithm are 4,5 and 6. split is applied on the associated dataflow relations → 4 , → 5 and → 6 . Each dataflow relation is split in three parts as depicted in Figure 2.(c). For → 5 : → 1 5 crosses hyperplane t (red), → 2 5 crosses hyperplane t + i (blue) and → 3 5 stays in a tile (green). This algorithm works pretty well for short uniform dependences → c : if fifo(c) before tiling, then, after tiling, the algorithm can split c in such a way that we get FIFOs. However, when dependences are longer, e.g. (t, i) → (t, i + 2), the target operations (t, i + 2) reproduce the tile execution pattern, which prevents to find a FIFO. The same happens when the tile hyperplanes are "too skewed", e.g. τ 1 = (1, 1), τ 2 = (2, 1), dependence (t -1, i -1) → (t, i). Figure 2.(d) depicts the volume of data to be stored on the FIFO produced for each depth. In particular, dotted line with k indicates iterations producing data to be kept in the FIFO at depth k. FIFO at depth 1 (dotted line with 1) must store N data at the same time. Similarly, FIFO at depth 2 stores at most b 1 data and FIFO at depth 3 stores at most b 2 data. Hence, on this example, each transformed channel requires b 1 + b 2 additional storage. In general the additional storage requirements are one order of magnitude smaller than the original FIFO size and stays reasonable in practice, as shown in the next section.

Experimental Evaluation

This section presents the experimental results obtained on the benchmarks of the polyhedral community. We demonstrate the capabilities of our algorithm at recovering FIFO communication patterns after loop tiling and we show how much additional storage is required.

Experimental Setup

We have run our algorithm on the kernels of PolyBench/C v3.2 [START_REF] Pouchet | Polybench: The polyhedral benchmark suite[END_REF]. Tables 1 and2 depicts the results obtained for each kernel. Each kernel is tiled to reduce I/O while exposing parallelism [START_REF] Alias | Data-aware Process Networks[END_REF] and translated to a PPN using our research compiler, Dcc (DPN C Compiler). Dcc actually produces a DPN (Data-aware Process Network), a PPN optimized for a specific tiled pattern. DPN features additional control processes and synchronization for Inria I/O and parallelism which have nothing with our optimization. So, we actually only consider the PPN part of our DPN. We have applied our algorithm to each channel to expose FIFO patterns. For each kernel, we compare the PPN obtained after tiling to the PPN processed by our algorithm.

Results Table 1 depicts the capabilities of our algorithm to find out FIFO patterns. For each kernel, we provide the channels characteristics on the original tiled PPN (Before Partitioning) and after applying our algorithm (After Partitioning). We give the total number of channels (#channel), the FIFO found among these channels (#fifo), the number of channels which were successfully turned to FIFO thanks to our algorithm (#fifo-split), the ratios #fifo/#channel (%fifo) and #fifo-split/#channel (%fifo-split), the cumulated size of the FIFO found (fifo-size) and the cumulated size of the channels found, including FIFO (total-size). On every kernel, our algorithm succeeds to expose more FIFO patterns (%fifo vs %fifo-split). On a significant number of kernels (11 among 15), we even succeed to turn all the compute channels to FIFO. On the remaining kernels, we succeed to recover all the FIFO communication patterns disabled by the tiling. Even though our method is not complete, as discussed in section 3, it happens that all the kernels fulfill the conditions expected by our algorithm (short dependence, tiling hyperplanes not too skewed).

Kernel

Before Partitioning After Partitioning #channel #fifo #fifo-split %fifo %fifo-split fifo-size total-size#channel #fifo fifo-size total-size trmm 2 Table 2 depicts the additional storage required after splitting channels. For each kernel, we compare the cumulative size of channels split and successfully turn to a FIFO (size-fifo-fail) to the cumulative size of the FIFOs generated by the splitting (size-fifo-split). The size unit is a datum e.g. 4 bytes if a datum is a 32 bits float. We also quantify the additional storage required by split channels compared to the original channel (∆ := size-fifo-split -size-fifo-fail / size-fifo-fail). It turns out that the FIFO generated by splitting use mostly the same data volume than the original channels. Additional resources are due to our sizing heuristic [START_REF] Alias | Bee+Cl@k: An implementation of lattice-based array contraction in the source-to-source translator Rose[END_REF], which rounds channel size to a power of 2. Surprisingly, splitting can sometimes help the sizing heuristic to find out a smaller size (kernel gemm), and then reducing the storage requirements. Indeed, splitting decompose channel into channels of a smaller dimension, for which our sizing heuristic is more precise. In a way, our algorithm allows to finds out a nice piecewise allocation function whose

Conclusion

In this paper, we have proposed an algorithm to reorganize the channels of a polyhedral process network to reveal more FIFO communication patterns. Specifically, our algorithm operates producer/consumer processes whose iteration domain has been partitioned by a loop tiling.

Experimental results shows that our algorithm allows to recover most of the FIFO disabled by loop tiling with almost the same storage requirement. Our algorithm is sensible to the dependence size and the loop tiling chosen. In the future, we plan to design a reorganization algorithm provably complete, in the meaning that a FIFO channel will be recovered whatever the dependence size and the tiling used. We also observe that splitting channels can reduce the storage requirements in some cases. We plan to investigate how such cases can be revealed automatically.

Figure 1 :

 1 Figure 1: Motivating example: Jacobi-1D kernel

Figure 2 :

 2 Figure 2: Impact of loop tiling on communication pattern Consider Figure 2.(c).Dependence 5 is partitioned in 3 parts: red dependences crossing tiling hyperplane φ 1 (direction t), blue dependences crossing tiling hyperplane t + i (direction t + i) and green dependences inside a tile. Since the execution order in a tile is the same than the original execution order (actually a subset of the original execution order), green dependences will clearly verify the FIFO communication pattern. As concerns blue and red dependences, source and target are executed in the same order because the execution order is the same for each tile and dependence 5 happens to be short enough. In practice, this partitioning is effective to reveal FIFO channels. In the next section, we propose an algorithm to find such a partitioning.

1 6 for each channel c 7 {→ 1 c 9 remove

 6719 split(→c,θP ,θC ) 2 for k := 1 to n 3 add(→c ∩{(x, y), θP (x) k θC (y)}); 4 add(→c ∩{(x, y), θP (x) ≈ n θC (y)}); 5 fifoize((P, C)) , . . . , → n+1 c } := split(→c,Pc,Cc); 8 if fifo(→ k c , ≺ θ Pc , ≺ θ Cc ) ∀k

Figure 3 :

 3 Figure 3: Our algorithm for partitioning channels

2 . 2

 22 Model at a Glance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Scheduling and Loop Tiling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 Polyhedral Process Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Table 1 :

 1 Detailed results

			1	2 50%	100%	256	512	3	3	513	513
	gemm	2	1	2 50%	100%	16	528	3	3	304	304
	syrk	2	1	2 50%	100%	1	8193	3	3	8194	8194
	symm	6	3	6 50%	100%	18	818	7	7	819	819
	gemver	6	3	5 50%	83%	4113	4161	7	6	4146	4162
	gesummv	6	6	6 100%	100%	96	96	6	6	96	96
	syr2k	2	1	2 50%	100%	1	8193	3	3	8194	8194
	lu	8	0	3 0%	37%	0	1088	11	6	531	1091
	cholesky	9	3	6 33%	66%	513	1074	11	8	788	1076
	atax	5	3	4 60%	80%	48	65	5	4	49	65
	doitgen	3	2	3 66%	100%	8192	12288	4	4 12289	12289
	jacobi-2d	10	0	10 0%	100%	0	8320	18	18	8832	8832
	seidel-2d	9	0	9 0%	100%	0	49952	16	16 52065	52065
	jacobi-1d	6	1	6 16%	100%	1	1153	10	10	1175	1175
	heat-3d	20	0	20 0%	100%	0 148608	38	38 158992 158992

Table 2 :

 2 Aliasfootprint is smaller than a single piece allocation. We plan to exploit this nice side effect in the future. Impact of splitting on storage requirements

	kernel	Size-fifo-fail Size-fifo-split	∆
	trmm	256	257	0%
	gemm	512	288 -44%
	syrk	8192	8193	0%
	symm	800	801	0%
	gemver	32	33	3%
	gesummv	0	0
	syr2k	8192	8193	0%
	lu	528	531	1%
	cholesky	273	275	1%
	atax	1	1	0%
	doitgen	4096	4097	0%
	jacobi-2d	8320	8832	6%
	seidel-2d	49952	52065	4%
	jacobi-1d	1152	1174	2%
	heat-3d	148608	158992	7%

RR n°9131


Mathieu Muniglia
email: mathieu.muniglia@cea.fr

Jean-Charles Le Pallec

Jean-Michel Do

Sébastien Verel

Design of a load-following management for a PWR reactor using an optimization method

In this paper, an optimization of the control rods parameters of a nuclear power plant is presented. Using a multi-physics simulator, the criteria of interest are computed in few minutes of computation and the optimization problem is solved thanks to a parallel asynchronous master-worker (1 + λ)-Evolutionary Algorithm. Two criteria are minimized independently: the first one based on the control diagram is a measure of the energy heterogeneity of the core, and a second one corresponding to the volume of effluents resulting from the use of the soluble boron. The main goal of this paper is then to analyze the solutions found by the optimization. Moreover, it is known that the burnup of the core has a strong influence on the use of the soluble boron, so different starting points have been considered, showing that the evolution of the performances of the good solutions as regards the two criteria are not the same. Finally, this paper shows that two gathered mono-objective optimizations are not satisfactory enough if one considers the burnup evolution, and that a bi-objective optimization must be considered to solve this problem.

I. INTRODUCTION

In the actual context of energetic transition, the increase of the intermittent renewable energies contribution (as wind farms or solar energy) is a major issue. On the one hand, the French government aims at increasing their part up to 30% [START_REF] Dumont | Ademe Energie 2030 : Production d'énergies renouvelables[END_REF] by 2030, against 6% today. On the other hand, their intermittent production may lead to an important imbalance between production and consumption. Consequently, the other ways of production must adapt to those variations, especially nuclear energy which is the most important in France. The power variations occur at different time scales (hour, day, or even week) and in order to counterbalance their effects on the electric grid, the nuclear power plants (NPP) are already able to adjust their production. NPPs which take part in the response of the power variations operate in the so-called load-following mode. In this operating mode, the power plant is controlled using control rods (neutron absorber) or soluble boron. However, the control rods may introduce unacceptable spatial perturbations in the core, especially if the power variations are large and/or fast, and the use of the soluble boron produces waste effluents that need to be processed. The purpose of this work is to optimize the manageability of the power plants to cope with a large introduction of intermittent renewable energies, and its final goal is to tune the control parameters in order to be able to make the load-following at a shorter time scale and larger power amplitude scale, or increase the safety margins to do so.

Such an optimization is a real challenge, considering the size of the search domain, the computation cost and the unknown properties of the objective function to minimize. Due to the design of the nuclear power plants, and in a goal to propose only simple modifications of the current management, 11 integer variables are used to describe the control rods such as speed, overlaps between rods, etc. (details are given in Sect. II.). The soluble boron use will be a consequence of those parameters selection. Therefore, the optimization is a large size combinatorial problem where no full enumeration is possible. Moreover, a multi-physics simulator is used to compute according to the variables of the main parameters of the control rods several criteria such as the evolution of the axial power offset, the rejected volume of effluent, etc.

The use of Evolutionary Algorithms (EA) in order to optimize some variables of a nuclear power plant as regards performance or safety is not new. Offline optimizations can already be found, and studies such as [START_REF] Arnaud | Selection Combinatory Algorithm for Loading Pattern Design of Light Water Reactor with Two Levels of Heterogeneity[END_REF] or [START_REF] De Moura Meneses | A new approach for heuristics-guided search in the In-Core Fuel Management Optimization[END_REF] deal with the In-Core Fuel Management Optimization (ICFMO) and loading pattern optimization which is a well-known problem of Nuclear Engineering and aims for instance at maximizing the use of the fuel (increase the cycle length for example) while keeping the core safe (minimize the power peak). Pereira and Lapa consider in [START_REF] Pereira | Coarse-grained parallel genetic algorithm applied to a nuclear reactor core design optimization problem[END_REF] an optimization problem that consists in adjusting several reactor cell variables, such as dimensions, enrichment and materials, in order to minimize the average peak-factor in a reactor core, considering some safety restrictions. This is extended in [START_REF] Sacco | Two stochastic optimization algorithms applied to nuclear reactor core design[END_REF] to stochastic optimization algorithms conceptually similar to Simulated Annealing. Sacco et al. even perform in [START_REF] Sacco | A Metropolis Algorithm applied to a Nuclear Power Plant Auxiliary Feedwater System surveillance tests policy optimization[END_REF] an optimization of the surveillance tests policy on a part of the secondary system of a Nuclear Power Plant, using a metaheuristic algorithm, which goal is to maximize the system average availability for a given period of time.

To the best authors' knowledge, the only optimizations of the plant operation are made online, like in [START_REF] Na | Design of a Fuzzy Model Predictive Power Controller for Pressurized Water Reactors[END_REF], where Na et al. develop a fuzzy model predictive control (MPC) method to design an automatic controller for thermal power control in pressurized water reactors. The objectives are to minimize both the difference between the predicted reactor power and the desired one, and the variation of the control rod positions. A genetic algorithm is then used to optimize the fuzzy MPC. Kim et al. propose in [START_REF] Kim | Design of a model predictive load-following controller by discrete optimization of control rod speed for PWRs[END_REF] another MPC by applying a genetic algorithm, to optimize this time the discrete control rod speeds. The offline Nuclear Reactor Operation Optimization (NROO) proposed in this paper and using an Evolutionary Algorithm (EA) is then an original combinatorial black-box problem. Moreover, the black-box objective function, which is also called fitness function in the context of evolutionary computation, is computationally expensive (one evaluation typically takes on average about 40 minutes) so that only few candidate solutions and their corresponding fitness value, which is the value of the criterion, can be computed. Otherwise, no property on the search space is a priori known.

With the increasing number of computing resources (cores, gpgpu, etc.), EAs which are population-based stochastic optimization algorithm, and more specifically parallel EAs become more and more popular to solve complex optimization problems. Usually, parallel EA [START_REF] Alba | Parallelism and evolutionary algorithms[END_REF] can be classified into two main models: the coarse-grained model (or island model) in which several EA share candidate solutions within the migration process, and the fine-grained model (or cellular model) where the search population is spread into a grid and evolutionary operators are locally executed. Moreover, the execution of parallel EA follows two main architectures. In the P2P architecture, the computation is fully distributed on all computation units, and in the Master-Worker (M/W) architecture, the fitness evaluation is on workers. The M/W architecture has been extensively used and studied [START_REF] Dubreuil | Analysis of a Master-Slave Architecture for Distributed Evolutionary Computations[END_REF]. It is simple to implement, and does not require sophisticated parallel techniques and communications. Two communication modes are usually considered. In the synchronous mode, the parallel algorithm is organized by round. The master sends candidate solutions on each worker for evaluation, and waits until receiving a response from all workers before the next round. In the asynchronous mode, the master does not need to wait, and communicates individually with each worker on-the-fly. The asynchronous mode could improve the parallel efficiency when the evaluation time of the fitness function vary substantially [START_REF] Wessing | Comparing Asynchronous and Synchronous Parallelization of the SMS-EMOA[END_REF]. In this paper, an asynchronous M/W parallel EA is used to solve the NROO problem.

The main goals of the paper are then : (i) perform for the first time an optimization of the control rods using an evolutionary algorithm (ii) analyze the efficient solutions as regards an axial offset criterion as well as the effluents volume and (iii) study the suitability of a bi-objective optimization.

The rest of this paper is organized as follows. The next section introduces the NROO problem, that is to say the tunable parameters and the criteria. The proposed algorithm is described in Sect. III., and in Sect. IV., the experimental analysis of the found solutions is conducted. At last, the paper concludes on the main results, and future works.

II. PROBLEM DEFINITION

The optimization process is based on the current loadfollowing transient [START_REF] Lokhov | Technical and economic aspect of load following with nuclear power plants[END_REF] and this analysis focuses on a single Pressurized Water Reactor (PWR) type (1300 MW) of the French nuclear fleet. The typical load-following transient considered here is presented on Fig. 1, and corresponds to the most penalizing scenario in the specifications of a nuclear power plant.

In a general way, when an electrical power variation occurs (demand of the grid) a chain of feedback is setting up in the whole reactor, leading to a new steady state. It is usual to take advantage of this self-regulation in the case of small variations, but the regulated variables such as the temperature or the pressure in the primary or secondary circuits may reach unacceptable values in case of such a load-following, possibly leading to damages of the whole system. The control rods are then used in order to cope with this variation, and maintain the primary coolant temperature close to the target. However, those control rods have to be handled carefully as they could cause axial or radial heterogeneity in the core, inducing high power peaks or Xenon oscillations.

Description of the System

The reactor core is a grid of square assemblies in a cylindrical vessel. Each assembly is 21 centimeters side length and about four meters high. There are 193 assemblies, split into two kinds : 120 assemblies made of Uranium oxide (UOX) and 73 ones made of Uranium plus Gadolinium oxides (UGd). Fig. 2 shows the core of the reactor : the UOX assemblies are in orange, the UGd ones in purple, and the green stands for the steel vessel. Each control rod is made of pins of a neutron absorber that are inserted together from the top of the core in some assemblies. The positions of the assemblies where they are inserted and the materials of which they are made correspond to the French "G" mode [START_REF] Lokhov | Technical and economic aspect of load following with nuclear power plants[END_REF]. There are two kinds of rods in this "G" mode : the black rods made of very absorbing pins (B 4 C and Ag+In+Cd) and the gray rods (giving their name to the mode), made of less absorbing pins (stainless steel and Ag+In+Cd). Those rods are then organized in two families, depending on their function : (i) the power shimming rods (PSR) and (ii) the regulation rods (TRR).

The first ones are used to shim the power effects during the power transient, and are split in four groups (four gray rods G1, eight gray rods G2, eight black rods N1 and eight black rods N2). All the rods of a same group move together, and the groups are inserted successively in this order : G1, G2, N1, N2, as it is shown in Fig. 3. An overlap is also defined between all the groups, so that they follow an insertion program as illustrated from frames (a) to (d). A convenient definition is used to describe the positions of the four groups using a single value : the totalizer (T). It represents the cumulative number of inserted steps, only incremented by the last moving group. Moreover, the totalizer is linked to the electrical power by a calibration function, established at the beginning of the cycle and updated every 60 equivalent full-power days to take into account the burnup effects. Those rods, however, do not enable to shim the evolution of the Xenon concentration during the transient, as well as the use of the soluble boron. The totalizer value (T ) is given on each frame, and the last moving group is in purple.

The second family is made of nine black rods gathered in a single group and enables a control of the average coolant temperature of the core (the targeted temperature, called reference temperature, is a linear function of the thermal power). Indeed, the Xenon evolution is managed mainly with the soluble boron, and another family of rods is then needed to cope with the resulting temperature discrepancies. This group moves independently and automatically, following a speed program depending on the difference between the reference temperature (T re f ) and the mean temperature (T m ) as shown in Fig. 4. One can see a dead band of ±0.8°C in which the rods do not move, avoiding continuous displacement and corresponding to the self-regulation of the core. Finally, as they are very efficient and for safety reasons, they can only move inside a maneuvering band of about 50 centimeters in the upper part of the core. For more details, please refer to [START_REF] Grard | Physique, fonctionnement et sûreté des REP[END_REF]. The dead band corresponds to the null speed and the maximal and minimal speeds (±72 steps/min) are for an absolute temperature difference larger than 2.8°C.

The variables to be tuned for the optimization are then the four nominal speeds and the three overlaps for the PSR, the maximal and minimal speeds, the dead band width and the maneuvering band height for the TRR. Eleven variables are then considered, and they are coded as integer values corresponding to a discrete number of steps or of temperature (the dead band is discretized by steps of 0.1°C). Tab. I summarizes the variables, their initial values (current management) and ranges. The values in the table take into account some technological and logical constraints. For example, the overlaps cannot be greater than the total height of the rods, the velocity ranges are bounded by the mechanisms, etc. A number of other variables could have been studied, like swaps between groups, or splitting groups, but the study is confined to the variables listed above for two reasons : simplify the problem for a first optimization, and be able to propose a solution without major technological breakthroughs and similar to the current one. Nevertheless, the search domain is huge (at least 3 × 10 20 possible configurations).

Criteria Definition and Calculation

This seek of simplification is even more understandable when it is known that the black-box evaluation function is very costly. Each unitarian calculation corresponds to a given management configuration running on a complete typical loadfollowing transient, corresponding to about eleven hours. The value of interest is then determined thanks to a model of the whole reactor described in [START_REF] Muniglia | A multi-physics PWR model for the load following[END_REF], and developed within the APOLLO3® [START_REF] Schneider | Apollo3 : CEA/DEN deterministic multipurpose code for reactor physics analysis[END_REF] calculation code. The optimization aims at minimizing this value of interest, which represents either a global operating criterion, based on the control diagram, or the rejected volume of effluents. The control diagram is used by the operator to manage the power plant and represents the evolution of the relative thermal power (P r ) as a function of the power axial imbalance ∆I:

∆I = P r × AO
where AO is the axial offset defined as AO = P T -P B P T +P B and standing for the unbalance between the lower and upper half parts of the core as regards the power. P T (resp. P B ) is the power in the upper (resp. lower) part of the core.

An example of such a diagram is to be found on Fig. 5, which draws the path of the state of the core during a power variation (blue line) and the bounds for this path. On the right side, the forbidden region (red line) is based on many studies and ensures the safety of the core in case of accidental situations. The impossible working region just comes from the definition (AO ∈ [-1, 1]). Finally, the green line starting at the same point as the path corresponds to a constant axial offset, and is called reference line in the following. The criterion derived from the control diagram is defined by:

C AO = 1 4 i |P 2 r,i+1 -P 2 r,i | • D(∆I i+1 ) + D(∆I i ) (1)
where

D(∆I i ) = |∆I i -∆I re f i |.
The point (P r,i , ∆I i ) represents the state of the core at the time step i in the control diagram, and ∆I re f i the power axial imbalance given by the reference line at the power P r,i . The criterion corresponds to the sum of all the areas as illustrated on Fig. 5, weighted by the relative power to take into account the fact that an important axial offset at high power is worse than at low power. Minimizing this criterion enables to reduce the area of the path and avoids being close to the forbidden region while staying close to the reference line.

The second criterion, V effluents , corresponds to the effluents volume produced during the transient. To increase (resp. decrease) the boron concentration in the core, some water of the primary circuit is replaced by highly boron-concentrated (resp. clear) water. Such a boration (resp. dilution) produces a lot of fluid wastes to be processed, called effluents. This volume is determined using the flow rate at which the new fluid is injected (equal to the one at which the current fluid is removed) by the pump :

V effluents = i Q i * ∆t i
where ∆t i and Q i are the length of the time step and the flow rate (in m 3 .s -1 ) of the pump at the time step i. The flow rate is given by :

Q i = - V ∆t i * log 1 + [B] i -[B] i-1 [B] i-1 -[B] * (2)
with V the volume of the primary circuit and [B] i the boron concentration at the time step i. [B] * is the concentration of the highly boron-concentrated water (7700 ppm) in case of boration, or of the clear water (0 ppm) in case of dilution.

Both criteria have to be minimized : one wants to reduce as well the axial perturbations on the core and the use of the soluble boron. They are a priori contradictory, because an important use of the control rods increases the axial perturbations, but reduces the use of the soluble boron. On the contrary, a limited use of the control rods produces few perturbations, but in that case, the soluble boron is used instead. The optimization process aims at finding a good compromise.

III. ASYNCHRONOUS PARALLEL EA

In this section, we propose a parallel EA dedicated to the NROO problem with computational expensive simulation cost, and to the high performance computing environment composed by a large number of computation units. The parameters of the EA algorithm are also studied.

Algorithm Design

Two keys guide the design of the optimization algorithm. First, the time of fitness evaluation computed by the simulator is on average about 40 minutes and with a large variance.

Second, a large number of computing units (w = 3072) are available to run the optimization algorithm, but they are available for only few hours (around 15 hours per job [16]). As a consequence, a master-worker (M/W) framework with fitness evaluation on each worker seems to be relevant. Indeed, on average we can expect one fitness evaluation every 0.78 second, meaning that the master node is not to be overflowed by the request of the workers, and with respect to the fitness evaluation time, an idle working time of few seconds will not reduce the performance. In addition, some simulations crash before the end of the calculation, increasing even more the discrepancies in calculation times. All considered, the model of the M/W has been made asynchronous : the workers are updated on the fly without a synchronization barrier, and each worker only computes the fitness value using the multi-physics simulator.

The number of evaluations per worker is small, on average 23 fitness function evaluations is possible on each worker within 15 hours of computation. As a consequence, the EA should converge quickly. We propose then an asynchronous (1 + λ)-EA where λ is the large number of computation units minus one. Roughly speaking, an iteration of the algorithm produces λ new candidate solutions (which are random variants of the so far best-known solution) by the mutation operator, and updates the best solution for the next iteration. The Algo. 1 and 2 show the details of the algorithm.

f i ← Msg[0] 9 if f i f * then 10 x ← x i ; f ← f i 11 end 12 x i ← Mutate(x ) 13
Send (non-blocking) Msg(x i ) to worker i

14 end 15 return x Algorithm 2: Asynchronous M/W (1 + λ)-EA on work- ers 1 Receive (blocking) Msg from master 2 x i ← Msg[0] 3 f i ← Evaluate x i with simulator 4 Send (blocking) Msg( f i ) to master
First, the algorithm on master node produces λ = w -1 quasi-random solutions (integer vectors of dimension n = 11) using a Design of Experiments (DoE) based on Sobol of quasirandom numbers. This initialization is used to improve the spreading of the initial solutions in the search space. Every initial solution is then sent asynchronously to a worker, running the Algo. 2 : receive the solution from the master, compute the fitness value by running the multi-physics simulator, and send back the result to the master node. In the meantime, the main loop of the Algo. 1 is executed on the master node : wait for a message from a worker i, and when the fitness value is received, the best so far solution is updated if necessary. Notice that the best solution is replaced by the new solution evaluated by the worker even when the fitness values are equals. In that way, the algorithm is able to drift on plateaus of the search space which are solutions with the same criterion value. A new candidate solution is then computed by the mutation operator of the best-known solution and sent in non-blocking mode to the same worker i. The master is then able to manage the requests of the other worker nodes by the asynchronous communication mode. The stopping criterium of the algorithm is the time limit.

The mutation operator is based on the classical mutation operator for vectors of numbers. The mutation rate p defines the parameter of the Bernouilli distribution to modify each number of the vector. Therefore, the number of modified variables follows a binomial distribution of parameters n and p, and the expectation of the number of modified variables is np. When an integer variable is modified according to the mutation rate, a random integer number is drawn using a uniform distribution centered on the current value. Let x j be the current value of the variable j, and δ j the gap defined by r.(ub jlb j ) where lb j and ub j are respectively the lower bound and the upper bound of the variable j defined in the Tab. I, and r ∈ [0, 1] is a mutation parameter. The new integer value of variable j after mutation is selected uniformly in the interval [x j -δ j , x j + δ j ] ∩ [lb j , ub j ] \ {x j }. The parameter r tunes the range width for the new value of variable after the mutation, and is expressed relatively to the total range width of the parameters.

In addition, to avoid multiple costly evaluations of the same candidate solution, a hash-map is used on the master node to save all evaluated solutions. The mutation is applied on the solution until a new candidate solution which is not in the hash-map is produced by the mutation random process.

Parameters Tuning

First, the performance of the algorithm with a baseline parameters setting is studied with 3072 computing units during 24 hours (approx. 73, 728 hours of CPU time). Following the value of the mutation rate parameter of 1/n commonly used in EA, the mutation rate has been set roughly to the inverse of the number of variables (p 0 = 0.1), so that the mutation operator modify on average one variable. The width of the random variation range has been arbitrarily set to about r 0 = 0.05 (5% of the total variation range of the variable).

The use of an asynchronous algorithm to avoid idle time is justified by the discrepancies of the computation costs from a candidate solution to another one : the mean computation time is 2426 seconds, and the faster computation is done in 1629 seconds whereas the longer is performed in 6169 seconds.

Fig. 6 shows the dynamic of the run. The C AO value is drawn as a function of the number of evaluations received by the master node. A point is plotted when the best solution so far is updated (included for equal fitness values) and the fitness values are normalized by the fitness value of the current management (see Tab. I). The solutions for which the number of evaluations is lower than 3072 are from the initial quasirandom population. Even if the baseline settings enable to reduce the fitness of about 40% compared to the current management, it can be seen that the number of strictly improving solutions is low (about 10 improving steps). The dynamic is a punctuated equilibrium dynamic with a lot of plateaus (made of solutions of the same quality), and therefore few improving solutions. For instance, the process is stuck on a plateau at the end of the run and almost 50, 000 fitness evaluations are necessary to find another plateau with strictly better solutions. This first experiment shows the relevance of the algorithm to find better solutions than the current management, but it also suggests that the settings of the mutation parameters could be improved.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.3 Four values of mutation rates p and mutation ranges r are then investigated : p ∈ {0.1, 0.2, 0.3, 0.4} and r ∈ {0.05, 0.1, 0.2, 0.5}. All the combinations are considered, given 16 possible mutation settings of the mutation operator. Moreover, to reduce the intrinsic random effect of the algorithm, each couple of mutation parameters values (p,r) have been launched five times with different initial populations generated by the Sobol sequence of quasi-random numbers. However, the five initial populations are the same for each couple of parameters settings. The algorithm is launched this time on 3072 computing units during five hours (approx. 15, 360 hours of CPU time per run), giving a total computation cost of more than 1, 2 × 10 6 hours of computation, and it was not possible to execute more than five runs per mutation settings.

Two sets of mutation parameters are shown on Fig. 6 (mut1 and mut2). The first set corresponds to (2p 0 , 2r 0 ), and is shown starting either from the same population as the baseline settings (init1), or from another pseudo-random population (init2). The second set is the best set found with this parametric study and corresponds to (3p 0 , 10r 0 ). The (mut1, init1) configuration is better than (baseline, init1), and enables to find another local minimum. In the same way (mut2, init2) is better than (mut1, init2). By increasing the mutation rate p and the random range r, one increases the differences between a solution and his mutant, improving the exploration phase. On the other hand, lower mutation rate and random range improve the exploiting around a good solution. Given the computing time, the number of cores, and the algorithm structure, it appears that the better the exploration, the better the solution. The initial population has also a huge influence on the solution, as it can be seen between (mut1, init1) and (mut1, init2).

One can see the importance of optimizing the parameters of the algorithm : a gain of about 65% is reached against 40% with the baseline settings.

IV. EXPERIMENTAL ANALYSIS 1. Mono-objective Optimizations

This section deeply analyzes the optimizations with respect to both criteria. To do so, the optimization process is launched again in order to optimize this time the effluents volume. The parameters of the algorithm are the one corresponding to the so far best configuration, which is (mut2, init2). The dynamic of the optimization is similar to the one for the axial offset criterion, and is not presented here. This optimization results and the ones coming from the best run as regards the axial offset criterion are then gathered (all the solutions are kept, and not only the improving ones as in Fig. 6), and the Pareto set is built. The Pareto set is made of all the non-dominated (or Pareto-optimal) solutions. Let X be the search space, and f : X → R n the fitness function which associates to each candidate solution the scalar values to be minimized (here n = 2). A solution x in search space is said dominated by a solution x if : ∀i ∈ {1, . . . , n} | z i ≤ z i and ∃ j ∈ {1, . . . , n} | z j < z j where z = f (x ) and z = f (x). A solution x is then nondominated if there does not exist any solution x that dominates it. The set of Pareto-optimal solutions is the Pareto set, and its mapping in objective space is the Pareto front, which is shown in Fig. ?? for the gathered optimization results.

It is satisfactory to see that the two mono-objective optimizations give a well-distributed front, and that the two optimizations contribute almost the same way to the front (same number of points from the optimization as regards C AO than V effluents ). The two criteria are contradictory as expected, since an optimization with respect to one degrades the other. Moreover, the current management is well improved for both criteria, and is dominated by about 15 points. It is interesting, however, to note that the Pareto front is centered around the current management, meaning that this management can be seen as an already improved solution as regards those two criteria.

Three candidate solutions are pointed out on this Pareto front : the one giving the best axial offset criterion (best AO), the one giving the best effluents volume (best Veff ), and one in the middle of the Pareto front (at (0.75,0.8)) called best q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.6 0.8 Compromise, giving one of the best compromise between the two criteria, and dominating the current management. Fig. 8 presents the corresponding control diagrams, and the evolution of the soluble boron concentration all along the load-following transient. It was said previously that the volume of effluents characterizes the use of the soluble boron, and consequently, the less the amplitude of the soluble boron concentration, the less the rejected volume of effluents.

First of all, Fig. 8 (top) shows that the area of the best AO path is well reduced (65% reduction) compared to the current management, and the axial offset stays very close to the reference line. In comparison, even if the area of the path of best Veff is small, the curve stays far from the reference line, explaining why C AO is important for this solution. As expected, best Compromise is between the two : reducing the use of the soluble boron but still ensuring a good behavior of the control diagram. The "S" shape of all the control diagrams excepted the best AO is due to the differential worth discrepancies of the PSR. One can avoid this effect by modifying their overlaps so that the reactivity insertion from one group to another is done in a smoother way, but also by reducing the velocities in such a way that the rods do not reach their expected position (given by the calibration function), reducing their axial effects on the core. Their final position is reached once finished the power variation. However, this strategy implies an important use of the soluble boron, which is used to shim the power effects (instead of the PSR), as it can be seen on Fig. 8 (bottom) from point A to B where the boron concentration increases. On the contrary, an important use of the PSR to shim the power effects like for best Veff, leads to a reduced use of the soluble boron (variation of 120 ppm in the case of best AO compared to 60 ppm for best Veff ).

Tab. II summarizes the parameters of the three previous candidate solutions, and remind the one of the current management. The overlaps of the best AO management are slightly higher than the current management ones, implying a priori higher perturbation on the core (a maximal overlap between two groups means that they are inserted together, contrary to a null overlap meaning that the first one is completely inserted before the second one starts moving). This is avoided by setting the G1 velocity to its minimal value, so that the rods cannot reach their expected position. Consequently, during the power decrease for instance, they are not enough inserted, meaning fewer effects on the core. Plus, the maneuvering band of the TRR is thin so they cannot be used neither. The boron is then used instead, and its almost uniform effect on the core explains why best AO ensures a path very close to the reference line, meaning that the axial offset is kept almost constant.

The best Veff configuration, however, maximizes the use of the PSR : G1 is inserted alone and when it reaches the bottom of the core, the power is around 75%PN. From that point, G2 and N1 are inserted together, and N2 starts moving shortly after them. This produces a huge decrease of the axial offset, as it can be seen on Fig. 8 (top). Moreover, the maneuvering band of the TRR is wide, so that the use of the soluble boron in this case is reduced to the minimum. Finally, the best Compromise configuration appears to be very close to the best AO, and the main difference is for the PSR velocities : with higher velocities, the power shimming is better, enabling to reduce the use of soluble boron. TABLE II. Summary of the parameters of the previous candidate solutions : current reference management (r.), best AO (a.), best Veff (b.) and best Compromise (c.).

PSR Overlaps

PSR Velocities TRR V. In order to generalize the comments made on the parameters of the three candidate solutions of the Pareto front, other candidate solutions are considered. For both optimization process (as regards C AO and V effluents ), all the solutions under a chosen normalized value are taken into account. This limit is set to 0.6 in case of C AO , and 0.75 in case of V effluents . The two values are different because the variation ranges on the two criteria are different, and the choice is made so that only the best solutions are kept, without removing too many solutions. Moreover, the number of kept solutions in both cases is almost the same (about 1350 solutions). Fig. 9 shows the parameter values distribution for the solutions found during the C AO optimization (left) and the V effluents optimization (right). q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 o1 o2 o3 v1 v2 v3 v4 V v mb db Normalized value q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 In both cases, the parameters o3, and v4 seem to have almost no influence on the values of interest, as the corresponding standard deviations are high. This can be explained by the fact that N2 is only used in few management configura-tions. Particular attention must be paid to the values of o1, o2, v1 and mb, which are the values that present the main discrepancies between the two optimizations. On Fig. 9 (left), the first quartile of the overlaps is always higher than the current management. For o1, the current management corresponds even to the lowest value. Moreover, the values of v1 and mb are close to the minimum. Consequently, the previous discussion on the strategy of the best AO management can be generalized to the good solutions found as regards C AO : by decreasing the velocity of the rods and the maneuvering band, one imposes the use of the soluble boron as the main control mechanism, and thus enable to improve the axial offset criterion. Keeping high overlaps opens, nevertheless, the possibility to decrease the use of the soluble boron by increasing the velocities. In the same way, the discussion on the strategy of the best Veff management is verified here, as it can be seen that o1 is minimal and o2 is maximal excepted for few solutions. On the one hand, the overlaps must be high in order to have a good shimming effect, as it is the case for best Veff, and leading to a reduced used of the soluble boron. The TRR enable to control the resulting average temperature and axial offset, and a high maneuveuring band is thus needed. But on the other hand, a too strong impact of the PSR on the core may also require the use of the soluble boron, which would increase the effluents volume, and that is probably why the first overlap is kept low.

o 1 o 2 o 3 v 1 v 2 v 3 v 4 V v

Influence of the Burnup

The soluble boron is also known to be used to counterbalance the decrease of reactivity all along the cycle, and its concentration never stops decreasing during the cycle. Usually, the maximum burnup for the load-following is given by a minimum amount of soluble boron left in the core. After a given burnup (about 80% of the cycle length), the reactor is not able anymore to load-follow, because the remaining soluble boron is not enough to counterbalance the Xenon effect. Moreover, the principle of boration and dilution of the soluble boron (Eq. ( 2)) implies that at low concentration, the boration is fast and the dilution is low. Consequently, at high burnup, a fast dilution is not possible, and that could also limit the possibilities of a core to load-follow. One can deduce that the burnup has an important effect on the performances of the solutions, especially if the use of the boron differs from a solution to another. As it is difficult to change online the management during the cycle, a candidate solution must be good all along the cycle to be chosen. Here, the robustness of the solutions with respect to the burnup is studied, by replaying the load-following transient for the solutions in the Pareto set, but with different initial states of burnup. The control diagrams of the best AO, best Veff and current management configurations are shown on Fig. 10, and correspond to a loadfollowing starting at 80% of the cycle length. The best Veff configuration is better than either the current management and best AO for both criteria. Its control diagram shape is the same as in Fig. 8 (top) : by minimizing the use of the soluble boron and maximizing the use of the control rods, one enables robustness as regards the burnup. When starting with a core at high burnup, the managements that need more soluble boron are consequently penalized. This is what Fig. 10 shows. For the current management and best AO, the needed boration and dilution speeds cannot be reached, making it difficult to fully counterbalance the power effects during the power decrease, and the xenon effect during the lower plateau. This is the cause of the drift observed during the lower plateau, resulting in a high C AO value. Knowing the importance of the burnup on the values of interest for those three management configurations, the effect of the burnup on the whole Pareto front is now shown in Fig. 11. From top to bottom, the initial burnup increases from 20% to 50% and then 80% of the cycle length (CL), and its effects on the front can be seen thanks to the color scale. The circle points stand for the initial Pareto front, and the triangle for the same management configuration, starting from the corresponding increased burnup. The references are at the center of the dotted cross.

As expected, the effect of the burnup on the Pareto front is important. Surprisingly, the more the burnup, the fewer the effects on the whole front. Even if the huge majority of the solutions is dominated by the initial Pareto front, some dominate it, meaning, as we discussed before that the effects are not uniform from a solution to another. More precisely, the color scale shows that the solutions optimized as regards the effluents volume tend to be improved as the burnup increases, whereas the solutions optimized as regards the axial offset criterion tend to be degraded. Indeed, the V effluents -optimized solutions start from the bottom right corner and go to the bottom left corner, meaning that they dominate the current management, whereas the C AO -optimized solutions start from the top left corner and go to the top right corner, meaning that they are dominated by the current management, and thus by the V effluents -optimized solutions. The solutions in the middle of the front stay close, most of the time, from the current management. The behavior of the Pareto front as a function of the burnup is still not well understood, but the final state (at q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.6 0.8

Normalized effluents volume

Burn-Up q 0*CL 0.2*CL q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.6 0.8

Normalized effluents volume

Burn-Up q 0*CL 0.5*CL q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.6 0.8 80% of the cycle length) corresponds to the discussion above. The absolute axial offset criterion of the V effluents -optimized solutions do not change much, whereas the value of the current management increases, making the normalized axial offset criterion decrease, and the normalized effluents volume of those solutions stays lower the one of the current management. Concerning the C AO -optimized solutions, the normalized effluents volume stays almost constant also (or slightly increases), but most of all, the axial offset criterion is increased, and gets closer to the reference. Fig. 10 shows indeed that the best AO and current managements are very similar to each other. One can conclude that the V effluents -optimized solutions are more robust than the C AO -optimized solutions as regards the burnup. It is not possible, however, to simply choose those solutions, because Fig. 8 shows that the maximal corresponding axial offset may be important, which is not acceptable with respect to the safety, particularly linear power and xenon oscillations. Consequently, the gathering of the two mono-objective optimizations is not enough (because not robust enough as it has been shown), and a bi-objective optimization is needed to find a compromise between V effluents -optimized solutions and acceptable axial offset variations.

V. CONCLUSIONS

A study concerning the optimization of the management of a Pressurized Water Reactor during a load-following transient have been presented, and the first results demonstrate the possibility of an improvement of the current management as regards two criteria which are the effluents volume and a criterion based on the axial offset calculation. It has been possible to improve the later by 65% and the former by 40%. However, some work is still to be done in order to optimize them jointly, and find good solutions as regards the two criteria. Those solutions will have, moreover, to be efficient as far as possible in the cycle length. This is the next step of the work : develop a multiobjective algorithm, based on the mono-objective one, enabling to find those solutions. The optimization will be done considering the beginning of the cycle in a first time, as it was the case here, and increased burnups in a second time. By comparing the results of those two bi-objective optimization, one will be able to point out the common solutions, thus efficient for the several burnups.

Fig. 1 .

 1 Fig.1. Load-following scenario considered for this optimization exercise. It is made of a power decrease at 5%PN/min to reach the lower plateau at 30%PN (point B), then the power is kept constant during 6 hours (until point C) before coming back to nominal power (point D).

Fig. 2 .

 2 Fig. 2. Position of the types of assemblies in a typical PWR1300.

Fig. 3 .

 3 Fig. 3. Insertion sequence of the Power Shimming Rods (PSR).The totalizer value (T ) is given on each frame, and the last moving group is in purple.

Fig. 4 .

 4 Fig.[START_REF] Pereira | Coarse-grained parallel genetic algorithm applied to a nuclear reactor core design optimization problem[END_REF]. Speed program of the Temperature Regulation Rods (TRR). The dead band corresponds to the null speed and the maximal and minimal speeds (±72 steps/min) are for an absolute temperature difference larger than 2.8°C.

Fig. 5 .

 5 Fig. 5. A control diagram showing the principle of the criterion computation to be minimized. The relative thermal power is expressed in percent of the nominal power (%PN).

Algorithm 1 : 7 Receive

 17 Asynchronous M/W (1+λ)-EA on master 1 for i in Workers do 2 x i ← Initialization using quasi-random numbers 3 Send (non-blocking) Msg(x i ) to worker i 4 end 5 f ← maximal value 6 while pending message and time is not over do Msg from worker i 8

Fig. 6 .

 6 Fig.6. Evolution of the normalized fitness as a function of the number of evaluations for some mutation parameters settings, and starting populations. Notice the log scale

Fig. 7 .

 7 Fig. 7. Pareto front obtained with the mono-objective optimizations as regards both criteria. The points coming from the two optimizations are distinguished. The values are normalized with the ones of the current management, which is at the center of the dotted cross.

Fig. 8 .

 8 Fig. 8. Control diagram and evolution of the soluble boron concentration for the three candidate solutions of the Pareto front, and for the current management. Current reference management (r.), best AO (a.), best Veff (b.) and best Compromise (c.). The A, B, C and D points corresponding to the power transient are indicated on the control diagram of the current management.

Fig. 9 .

 9 Fig. 9. Summary of the parameter values for the best solutions found with the C AO (left) and V effluents (right) optimization. The bottom (resp. top) of the boxes corresponds to the first (resp. third) quartiles, and the lower (resp. upper) whiskers extend from the boxes to the lowest (resp. highest) value within 1.5 times the distance between the first and third quartile. Values beyond or above are plotted as points. The median values are represented by the horizontal segments inside the boxes. The red points correspond to the current management.

Fig. 10 .

 10 Fig. 10. Control diagrams for the current management, best AO and best Veff at 80% of the cycle length.

Fig. 11 .

 11 Fig. 11. Effects of an increasing burnup on the Pareto front. 0.2CL (top), 0.5CL (middle) and 0.8CL (bottom). The values corresponding to managements starting from an increased burnup are normalized with respect to the values of the current management, starting from the same increased burnup.

TABLE I

 I

		. Integer variables of the design : lower bound (l.),
	upper bound (u.), and values of the current reference manage-
	ment (r.). The dead band (db) parameter is expressed in tenth
	of degree, and all the other variables are expressed in steps.
	PSR Overlaps		PSR Velocities	TRR V. TRR B.
	o 1 o 2 o 3	v 1	v 2	v 3	v 4	V	v	mb db
	l. 0	0	0	10 10 10 10 42 3	7	8
	u. 255 255 255 110 110 110 110 102 13 117 16
	r. 70 80 95 60 60 60 60 72 8	27 8

VI. ACKNOWLEDGMENTS

APOLLO3 is a registered trademark of CEA. The authors would like to thank IDEX and PS2E for their financial support. The authors gratefully acknowledge AREVA and EDF for their long-term partnership and their support. This work was performed using HPC resources from GENCI-TGCC (Grant 2016-t2016107616).

The authors would especially like to thank Hubert Grard from INSTN for his help in the setting of the simulator as well as for his precious advises concerning the analysis of the results.


Alain Ghio

Parallel distributed

come