Yann Barsamian
email: ybarsamian@unistra.fr

Arthur Charguéraud
email: arthur.chargueraud@inria.fr

Sever A Hirstoaga
email: sever.hirstoaga@inria.fr

Michel Mehrenberger

Efficient Strict-Binning Particle-in-Cell Algorithm for Multi-Core SIMD Processors

Keywords: Particle-in-cell, Plasma physics, Multi-core, SIMD architecture, Shared memory, Chunks, Strict binning, Magnetized electron holes

Particle-in-Cell (PIC) codes are widely used for plasma simulations. On recent multi-core hardware, performance of these codes is often limited by memory bandwidth. We describe a multi-core PIC algorithm that achieves close-to-minimal number of memory transfers with the main memory, while at the same time exploiting SIMD instructions for numerical computations and exhibiting a high degree of OpenMPlevel parallelism. Our algorithm keeps particles sorted by cell at every time step, and represents particles from a same cell using a linked list of fixed-capacity arrays, called chunks. Chunks support either sequential or atomic insertions, the latter being used to handle fast-moving particles. To validate our code, called Pic-Vert, we consider a 3d electrostatic Landau-damping simulation as well as a 2d3v transverse instability of magnetized electron holes. Performance results on a 24-core Intel Skylake hardware confirm the effectiveness of our algorithm, in particular its high throughput and its ability to cope with fast moving particles.

Introduction

The Particle-in-Cell (PIC) method enables large-scale simulations of plasma physics. PIC simulations are, for example, key to the design of ITER fusion reactor, and they also apply to other domains, e.g., astrophysics. As of 2018, PIC simulations accommodate at the order of 10 13 particles, involving the hundreds of thousands of cores available on the world's top super-computers.

To reach such a scale, state-of-the-art PIC codes exploit parallelism available at three levels: inter-node parallelism using, e.g., MPI; shared-memory multi-threading using, e.g., OpenMP, and register-level parallelism using SIMD instructions. A number of implementations also leverage GPUs or MICs. Implementations include EMSES [START_REF] Nakashima | Large Scale Manycore-Aware PIC Simulation with Efficient Particle Binning[END_REF], GTC-P [START_REF] Tang | Extreme Scale Plasma Turbulence Simulations on Top Supercomputers Worldwide[END_REF], ORB5 [START_REF] Jocksch | A Bucket Sort Algorithm for the Particle-In-Cell Method on Manycore Architectures[END_REF], OSIRIS [START_REF] Fonseca | Exploiting multi-scale parallelism for large scale numerical modelling of laser wakefield accelerators[END_REF], PI-CADOR [START_REF] Surmin | Co-design of a Particle-in-Cell Plasma Simulation Code for Intel Xeon Phi: A First Look at Knights Landing[END_REF], PIConGPU [START_REF] Bussmann | Radiative Signatures of the Relativistic Kelvin-Helmholtz Instability[END_REF], PSC [START_REF] Germaschewski | The Plasma Simulation Code: A modern particlein-cell code with patch-based load-balancing[END_REF], VPIC [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF].

A recent paper [START_REF] Tang | Extreme Scale Plasma Turbulence Simulations on Top Supercomputers Worldwide[END_REF] studied GTC-P performances in details, and points out that: "metrics such as flop/s or percentage-of-peak are less relevant for the predominantly memory-bound gyrokinetic PIC methods." The authors then present a model able to predict execution time as a function of data transfers. Most predominant is the intra-node operations on the shared memory (60% to 80% of the execution time). Their cost is decomposed between in-cache accesses, contiguous accesses, and random accesses-the latter being the most costly. This study shows that, to improve the performance of multi-core (intra-node) processing in PIC simulations, we must decrease the amount of costly memory accesses. Of course, we must do so by preserving the OpenMP-level parallelism as well as the crucial use of SIMD instructions.

The strict-binning approach to implementing the PIC method enables significant reduction in the number of random accesses and cache misses [START_REF] Tskhakaya | Optimization of PIC codes by improved memory management[END_REF][START_REF] Durand | A Packed Memory Array to Keep Moving Particles Sorted[END_REF][START_REF] Nakashima | Large Scale Manycore-Aware PIC Simulation with Efficient Particle Binning[END_REF][START_REF] Barsamian | A Space and Bandwidth Efficient Multicore Algorithm for the Particle-in-Cell Method[END_REF]: at every time step, particles that fall in the same cell are stored together. Doing so brings two main benefits. First, the electric field values can be read only once per cell, avoiding numerous cache misses and allowing SIMD computations when updating particle velocities. Second, the representation saves the need to store a cell index for each particle, thereby saving memory loads and stores.

A central challenge with the strict-binning approach is the representation of the dynamically-sized bins storing the particles. A state-of-the-art proposal by Nakashima et al. [START_REF] Nakashima | Large Scale Manycore-Aware PIC Simulation with Efficient Particle Binning[END_REF] organizes particles in a big array, ordering them according to their cell index, and leaving variable-size gaps between the groups of particles associated with each cell. Yet, this approach suffers from two important limitations. First, as particles move, maintaining the variable-size gaps requires costly operations for shifting particles. Second, the algorithm, which uses a coloring scheme [START_REF] Kong | Particle-in-cell simulations with charge-conserving current deposition on graphic processing units[END_REF] to avoid data races when processing the cells in parallel, does not handle well fast-moving particles (particles moving more than a couple cells away at a given time step): it resorts to sequential processing for these particles.

These two limitations are exacerbated when the percentage of crossing particles (particles changing cells at each time step) increases, to the point of possibly becoming a major bottleneck. For example, in a parallel execution using 64 cores, having as little as 0.5% of fast-moving particles can result in a 32% slowdown on the total execution time due to the sequential processing of these particles. 1 We propose an algorithm implementing strict-binning for the PIC method that addresses the aforementioned limitations, while still supporting efficient OpenMP/SIMD parallelization of all critical loops. Our algorithm leverages the use of chunk bags, i.e. linked lists of fixed-capacity arrays, to achieve SIMDfriendly storage of particles with limited memory overheads. These chunk bags are furthermore devised to support atomic push operations, which are used to handle fast-moving particles within the main parallel loop. Our algorithm minimizes the amount of memory transfers: at each time step, each particle gets read from and written to memory exactly once. In particular, no further move or reordering is ever required, regardless of the percentage of fast-moving particles. This algorithm is efficient provided that the average number of particles per cell exceeds a couple hundreds. Although laser-driven particle acceleration simulations can use as few as 30 particles per cell [START_REF] Bussmann | Radiative Signatures of the Relativistic Kelvin-Helmholtz Instability[END_REF], large-scale, high-precision simulations may involve hundreds to thousands of particles per cell [START_REF] Tskhakaya | Optimization of PIC codes by improved memory management[END_REF][START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF][START_REF] Germaschewski | The Plasma Simulation Code: A modern particlein-cell code with patch-based load-balancing[END_REF].

Through the rest of this paper, we describe our algorithm, comment on its theoretical properties (space usage, parallelization of critical loops, amount of memory transfers), discuss performance results (bandwidth usage, impact of fast-moving particles), numerical results (simulation of Landau-damping and of transverse instability of magnetized electron holes), and related work.

2 An Efficient, Strict-Binning, Multicore PIC Algorithm Fig. 1 describes the PIC method, applied to the resolution of the Vlasov-Poisson system, which models the time evolution of the distribution function of charged particles in a plasma. Following the Cloud-in-Cell model [3, Sect. 2.6.], we interpolate the electric field and accumulate the charges linearly from/to the eight corners of the grid cell where each particle lies. The Poisson solver takes less than 5% of the execution time, we thus focus our attention on the particle loop.

Our implementation performs all computations in double precision, with the exception of positions, which are stored using the "index plus offset" representation [4, III.E.], whereby the position of a particle is described relative to the corner of the cell containing the particle, using 3 float values, yielding sufficient precision. In the strict binning approach, the index of the containing cell is implicit. Thus, each particle admits a 36-byte representation.

Particles are stored in fixed-capacity arrays (chunks, e.g., [START_REF] Barsamian | A Space and Bandwidth Efficient Multicore Algorithm for the Particle-in-Cell Method[END_REF]). Several chunks might be needed to store all the particles contained in a same cell. Each cell is thus described by a linked list of chunks (a chunk bag). The number of particles in a chunk, denoted by K thereafter, should enable efficient vectorization (K is a multiple of 16 for 512-bit registers), and at the same time be large enough to tame the cost of following a pointer from a chunk to the next (e.g., 128 or 256).

Benchmarking on the hardware considered reveals that the structure of arrays (SoA) layout, which enables better vectorization, improves performance compared with the array of structures (AoS) layout. The memory layout we use for the particles is summarized on the next page. (Arrays should be aligned.) The chunk bag data structure supports O(1) insertion of a particle (adding a fresh chunk if needed), O(1) merge of two bags thanks to the back field (note that chunk compaction is not needed), and O(n) iteration over the contents, all with excellent constant factors. Furthermore, unlike chunks introduced by prior work [START_REF] Barsamian | A Space and Bandwidth Efficient Multicore Algorithm for the Particle-in-Cell Method[END_REF], our chunk bags are devised to support a thread-safe atomic insertion operation. Atomic insertions are central to the handling of fast-moving particles. Atomic insertion uses a fetch-and-add instruction to reserve a slot in the chunk for the particle. If a thread attempts to reserve the one-past-the-end slot, it acquires responsibility to extend the bag with a fresh chunk, in which case it sets the next pointer of the fresh chunk to the current front pointer of the bag, and sets the front pointer of the bag to the address of the fresh chunk.

struct
When processing a chunk of particles, the algorithm first updates velocities and positions, then migrates the particles to different chunks, depending on the cell associated with their new position. Once all particles from the chunk are processed, the chunk is stored into a (per-core) free list, so as to be subsequently reused to extend a bag whose last chunk becomes full. Our algorithm preserves the following invariant: at the beginning of a time step, all the particles are stored in at most N/K + 2 • nbCells chunks, where N denotes the total number of particles, and K denotes the number of particles per chunk.

To dispatch particles according to their target cells, we associate two bags with each cell: a private bag, accessed at most by one thread at a time; and a shared bag, accessed concurrently, to handle fast-moving particles. To initialize these two bags, we need an additional 2 • nbCells empty chunks. In total, we need N/K + 4 • nbCells chunks. We have proved that this number of chunk suffices at any point of a simulation, regardless of how particles move. Thus, the space used by our algorithm, in addition to the minimal amount of memory needed to represent the particles, grows in proportion with 4K • nbCells • sizeof(particle). 2In order to maximize the number of insertions into private bags while preserving a high degree of OpenMP parallelism, we follow the coloring scheme proposed by Kong et al. [START_REF] Kong | Particle-in-cell simulations with charge-conserving current deposition on graphic processing units[END_REF], and generalized from 2d to 3d by Nakashima et al. [START_REF] Nakashima | Large Scale Manycore-Aware PIC Simulation with Efficient Particle Binning[END_REF]. The idea is to fill the space with tiles, of size 2x2x2 (or more), in a regular manner. Tiles are colored using 8 different colors in such a way that two adjacent tiles have distinct colors. At each of the 8 color phases, 1 8 of the tiles are processed, in parallel by nbCores threads. 3 Because cells processed in parallel by distinct threads are at least 2 cells away from each other, all the particles that bag particles[0..nbCells-1]; // Particles by cell, at current time step bag particlesNextPrivate[0.

.nbCells-1], particlesNextShared[0..nbCells-1]; double ρ[0..X][0..Y][0..Z], E[0..X][0..Y][0..Z]; double ρNext[0..nbCores-1][0..nbCells-1][0..7]; //
ρ[x][y][z] += ρNext[idCore][idCell][i] ρNext[idCore][idCell]
[i] = 0 Compute E from ρ using a Poisson solver and set ρ to 0 // FFTW + OpenMP Fig. 2. Our parallel algorithm for the PIC method on multicore architectures. move, at a given time step, no more than one cell away (no more than half a tile away, in general) can be pushed into private bags, in a thread-safe manner.

The pseudo-code of our algorithm appears in Fig. 2. Particles from a same cell are processed sequentially by a same thread. To benefit from SIMD performance, we split the loop over each chunk. (If one chunk does not fit into L1 cache, additional splitting is needed.) First, our code updates velocities (line 12). Second, it computes the new positions (line 14), introducing an auxiliary array for storing the new cell indices. Third, it sequentially pushes each particle into the chunk associated with its target cell. If the target cell lies in the current tile, or lies in the closer half of an immediate neighboring tile, a non-atomic insertion is performed on a private bag (line 19). Otherwise, an atomic insertion is performed on a shared bag (line 21). Note that the boolean condition involved can be evaluated using a simple arithmetic test.

Once all the particles are processed, the algorithm merges, for each cell, its private bag with its shared bag (line 26). No chunk compaction is performed at this point; as a result, the bag associated with one cell may contain up to 2 nonfull chunks. Thus, there are at most N/K + 2 • nbCells nonempty chunks at the beginning of the next time step. It follows that at least 2 • nbCells empty chunks must have been freed during the current time step. This number corresponds exactly to the number of chunks needed to initialize the private and the shared bags for the next time step. Our algorithm performs this initialization efficiently in parallel (using a prefix sum array, based on the sizes of the per-core free lists).

We next describe the treatment of the charge density and the electric field (ρ and E). When processing particles from one cell, the algorithm first reads from memory the values of the electric field on the 8 corners of that cell (line 9). Importantly, thanks to the strict-binning approach, this data needs only be loaded once from memory. As particles are processed and moved to their target cells, the charge of each particle is accumulated (line 22) into the array ρNext, which, at the end of the time step, is used to update E for the next iteration. We exploit a recently-proposed, ingenious technique allowing to accumulate the charge on the 8 corners using SIMD instructions [START_REF] Vincenti | An efficient and portable SIMD algorithm for charge/current deposition in Particle-In-Cell codes[END_REF]. Concretely, the array ρNext involves some amount of redundancy: for each cell, 8 values are stored adjacently in memory, describing the charge on the 8 corners of that cell. At the end of a time step, the charge at a grid point is computed by summing the values associated with the 8 cells that have this grid point as one of their corners (line 32).

We considered two different possibilities for updating ρNext. The first possibility is to decompose ρNext into a private array and a shared array, just like we do for bags of particles. In this approach, only the deposit of the charge of fast-moving particles triggers atomic operations; for all others particles, we can use SIMD operations. The second possibility is to decompose ρNext into nbCores arrays. In this approach, each core has exclusive access to its charge array, so all accesses use SIMD operations. The downside is a slight increase in the memory footprint, and in the time needed to sum up the values. However, under our assumption of a reasonably large number of particles per cell, these additional costs are tiny in front of the gains. Thus, we opted for the latter approach.

Under the assumption of (at least) hundreds of particles per cell in average, the operations for manipulating chunks (following pointers, pushing/popping in free lists) and for manipulating per-cell information are all well amortized. Overall, our algorithm is not far from optimal in terms of memory transfers.

Optimization when particles move at most one cell per time step. For simulations whose physical parameters ensure that movement is restricted to immediate neighboring cells (e.g., [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF][START_REF] Germaschewski | The Plasma Simulation Code: A modern particlein-cell code with patch-based load-balancing[END_REF]), we can optimize our algorithm by removing the shared bags altogether. In this case, we require only N/K + 2 • nbCells + nbCores chunks, and do not need any atomic insertion operation. Likewise, ρNext can be stored in a single array (indexed by cells and by corners).

To assess correctness and performance of our code, which is called Pic-Vert, we considered two classical test cases: a 3d Landau-damping simulation and a 2d3v electron hole simulation. Section 4 presents details on these experiments, and argues that the numerical results produced by our simulation match the expected results. In the remaining of this section, we discuss performance results.

Experimental hardware is an Intel Xeon Platinum 8160 @ 2.1 GHz (Skylake), with 96 GB of RAM, 6 memory channels, and 24 cores. Our C code was compiled using Intel C Compiler 17.0.4, and the FFTW3 library [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF] for the Poisson solver.

The algorithm depends on two parameters. First, we use tiles of size 2x2x2 for the coloring. Tiles of size 4x4x4 lead to similarly good performances. Using larger tiles degrades performance. Second, we use K = 256 for the chunk capacity. Larger values of K increase the space usage and do not reduce the execution time. Smaller values of K increase the execution time overheads: +12% for K = 128, and +52% for K = 64. Note that, for K = 256, the memory "slack", which is equal to 4 • nbCells • sizeof(chunk), represents in the Landau-damping simulation only 13% of the amount of memory strictly required for representing the particles.

Achieved throughput. For the end-user of a simulation, the metric that matters is the number of particles processed per second. The Pic-Vert code achieves:

-740 million particles per second (30.8m/s/core) in the 3d Landau-damping simulation, where 31% of the particles change cell at each iteration; -910 million particles per second (37.9m/s/core) in the 2d3v electron hole simulation, where 32% of the particles change cell at each iteration.

Analysis in the roofline performance model. As argued in Sect. 2, our algorithm performs not far from the minimal number of memory operationsa key feature for PIC simulations hit by the memory bandwidth bottleneck. With this property in mind, it is interesting to compare the memory bandwidth achieved by our algorithm against the capacity of the hardware. Consider the Landau-damping simulation. The memory bandwidth achieved is 53.6 GB/s. 4 The theoretical peak advertised by the manufacturer is 127.99 GB/s. The Stream benchmark [START_REF] Mccalpin | Memory Bandwidth and Machine Balance in Current High Performance Computers[END_REF], which aims at evaluating the practical peak using a few microbenchmark programs, and which is commonly used as a baseline, provides the measure 98.2 GB/s. Our algorithm thus achieves 42% of the theoretical peak and 55% of the practical peak bandwidth. Reaching higher percentage in a PIC simulation appears to be very challenging.

Our algorithm is memory bound. In general, an algorithm may be compute bound (i.e. limited by the number of floating-point operations per second) or 4 The bandwidth is obtained by multiplying the size of a particle (36 bytes, plus 64 K bytes to account for chunk headers) by the number of particle processed per second (740 million), and by a factor 2 (one read plus one write). It would be very interesting to compare with other algorithms. Unfortunately, such numbers are rarely advertised, and comparisons are often hazardous due to differences in hardware, in particle representation (e.g., float vs double), in numerical schemes, etc. memory bound (i.e. limited by the number of bytes per second transfered from main memory) depending on its operational intensity, defined as the number of operations performed divided by the number of bytes moved from or to the main memory. We computed the operational intensity of the 3d code by counting the number of floating point operations per particle (79 operations in single-precision and 65 in double-precision, which leads to 209 operations when normalized to single-precision), and counting the number of bytes used to represent a particle (36 bytes, plus 0.25 byte to account for chunk headers). We thus derive that our 3d code has an operational intensity equal to 209/(2 • 36.25) ≈ 2.9. Similarly, we computed the operation intensity for the 2d3v code to be 114/(2 • 32.25) ≈ 1.8. Fig. 3 represents the bounds on computation and memory bandwidth, in a chart showing the operational intensity on the x-axis, and the computation performance on the y-axis [START_REF] Williams | Roofline: An Insightful Visual Performance Model for Multicore Architectures[END_REF]. Note that both axes are log-scale. The computation bound is an horizontal line, at 1,612 GFlops/s (billion floating-point operations per second), a figure provided by the hardware manufacturer. The theoretical and practical memory bounds (bytes/s) are diagonal lines, because the bound in performance (flop/s) is equal to the operational intensity (flop/byte) multiplied by the memory bandwidth (bytes/s). Each diagonal line meets the horizontal line at the point of break-even between memory bound and compute bound.

Efficient processing of fast-moving particles. In addition to being memory efficient, our algorithm also benefits from another key feature not found in prior strict-binning algorithms: fast-moving particles are handled efficiently within the main parallel loop. For a particle moving more than half a tile away, we require only one extra atomic operation. Moreover, the contention associated with this atomic operation is relatively limited. Indeed, for two atomic operations to be issued on the same memory cell at the "same time" (i.e., close enough in time for a race on the cache line to occur), it must be the case that two particles taken Particles that move 1 cell away 8.0% 8.0% 8.0% 8.0% 8.0% 8.0% 8.0% Particles that move 2 cells away 0 0.7% 1.9% 3.1% 4.3% 5.6% 4.4% Particles that move 3 cells away 0 0 0 0 0 0.2% 1.4% Particles pushed atomically (line 21) 0.0% 0.4% 1.0% 1.6% 2.2% 3.1% 3.7% Slowdown w.r.t. first column 0 0.0% 0.9% 3.8% 4.4% 4.2% 7.0%

Table 1. Impact on performance of increasing the percentage of fast particles.

from two distinct tiles spaced away by at least one full tile are moving towards the same cell of a third tile, at the "same time". Thus, the performance of our algorithm should be relatively independent form the heat.

To empirically evaluate the impact of fast-moving particles, we consider a simulation in which we artificially varied the initial distribution of particle velocities. To that end, we manually tuned these distributions in such a way as to obtain several test cases with increasing number of fast-moving particles. 5 Each test case is reflected by a column from Table 1. More specifically, the three first rows show the percentage of particles that move away from 1, 2 or 3 cells from their current grid cell at each time step (no particle move further away).

By instrumenting the code, we measured the number of push operations that trigger an atomic write (line 21 from Fig. 2). These numbers, relative to the total number of particles, appear in the fourth line of the table: they vary from 0% to 3.7%. The last row of Table 1 gives the corresponding slowdown on the total execution time. Figures show that even when the percentage of particles whose move require an atomic operation is as high as 3.7%, the cost of processing these fast moving particles remains fairly limited: +7.0%. In comparison, any alternative algorithm that sequentially processes 3.7% of the particles in a 24core execution would suffer at least from a +85% slowdown (recall Sect. 1).

Scaling. Although inter-node parallelism is orthogonal to the focus of the present paper, we used particle decomposition to scale our algorithm on 128 Skylake sockets (each with 24 cores, 12.3 TB of RAM in total), using one MPI process per socket. We simulated Landau-damping with 256 billion particles, achieving a throughput of 89.6 billion particles per second: 123x speedup w.r.t. one socket.

Numerical Results

3d3v Landau-damping. We consider a classical Landau-damping test case [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF], simulating 2 billion particles on a 64×64×64 grid, for 500 time steps. We use the same parameters as in [START_REF] Ricketson | Sparse grid techniques for particle-in-cell schemes[END_REF]: time step of 0.05, periodic boundary conditions on spatial domain Ω = [0, 22] 3 and initial distribution function:

f 0 (x, y, z, v) = 1 (2π) 3/2 e -|v| 2 /2 L(x)L(y)L(z)
with L(w) = 1 + 0.05 cos(πw/11). Fig. 4 represents the evolution of electric energy. It shows that the decay slope in our simulation is in accordance with the theoretical value γ = -0.008466 obtained from the dispersion analysis. 2d3v Electron hole. We consider a more complex test case proposed by Muschietti et al. [START_REF] Muschietti | Transverse Instability of Magnetized Electron Holes[END_REF]. We simulate 64 billion particles on a 512 × 512 grid (on 32 Skylake sockets). Time step is 0.1 and spatial domain is [0, L] 2 , with L = 32. The initial function is f (x, y, v) = F 1 (v 2 1 -2φ(x, y))e -50(v 2 2 +v 2

3 ) with potential φ(x, y) = e -0.5((x-L/2)/∆ -0.3 cos(0.39y)) 2

, with ∆ = 3 and F 1 defined as:

F 1 (w) =    √ -w π∆ 2 1 + 2ln( ψ -2w ) + 6+( √ 2+ √ -w)(1-w) √ -w π( √ 2+ √ -w)(4-2w+w 2 ) , for -2ψ ≤ w < 0, 6 √ 2
π(8+w 3 ) , for w > 0.

The external magnetic field is aligned with x and has amplitude B 0 = 0.2. Fig. 5 shows the charge density ρ(t, x, y) = 1 -f (t, x, y, v)dv, on the left at time t = 20, and on the right at time t = 40. These results are qualitatively similar to those from Muschietti et al. [START_REF] Muschietti | Transverse Instability of Magnetized Electron Holes[END_REF].

In addition, we studied the convergence of the simulation with respect to the number of particles and to the grid size. To that end, we compare, for different settings of these two parameters, the time evolution of a quantity representative of the instability. 6 Results appear in Fig. 6. They show that using a small 32×32 grid with 200 million particles as considered by Muschietti et al. exhibits the correct qualitative behavior up to t = 50, but diverges beyond this point.

For a quick simulation, it appears preferable to use a 64×64 grid with only 20 million particles, as it gives quantitatively accurate results up to t = 50. For longer simulations, our results show that using a 128×128 or a 256×256 grid with 200 million particles suffices to give accurate results up to t = 100. Indeed, the two corresponding curves are close to that of our large-scale simulation, which uses a 512×512 grid with 64 billion particles (the top-most curve at t = 100).

Technical Comparison to Related Work

We organize the discussion of related work by focusing on three main criteria: strict or non-strict binning, representation of particles, and treatment of data races arising when two threads push data onto a same target cell.

To ensure efficient accesses to the electric field and charge arrays, locality is essential in PIC simulations. In numerous algorithms, particles are stored in an array and ordered by their cell index. Yet, because particles move in the grid, their locality in the array decreases at each time step. Thus, reordering operations must be performed: either at every time step, to maximize locality (e.g., [START_REF] Jocksch | A Bucket Sort Algorithm for the Particle-In-Cell Method on Manycore Architectures[END_REF]), or at some lower frequency, reducing locality but mitigating the cost of re-sorting (e.g., [START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF][START_REF] Germaschewski | The Plasma Simulation Code: A modern particlein-cell code with patch-based load-balancing[END_REF]). Depending on the option, performances suffer either from numerous costly random accesses or from suboptimal locality.

Rather than sorting, other algorithms rely on coarse-grained binning [START_REF] Bussmann | Radiative Signatures of the Relativistic Kelvin-Helmholtz Instability[END_REF][START_REF] Tang | Extreme Scale Plasma Turbulence Simulations on Top Supercomputers Worldwide[END_REF][START_REF] Surmin | Co-design of a Particle-in-Cell Plasma Simulation Code for Intel Xeon Phi: A First Look at Knights Landing[END_REF][START_REF] Vincenti | An efficient and portable SIMD algorithm for charge/current deposition in Particle-In-Cell codes[END_REF]. Particles are organized in super-cells (of size, e.g., 10x10x10), and a dynamically-sized data structure is used to represent particles from a same supercell. For example, attribute tiles [START_REF] Bussmann | Radiative Signatures of the Relativistic Kelvin-Helmholtz Instability[END_REF] have been used to store particles on GPU using doubly linked lists of fixed-capacity arrays. Unlike with our chunks, particles in an attribute tile are processed in place. If a particle moves to a different super-cell, it is migrated to a transfer buffer, and its slot is marked as a hole in an auxiliary bitmap. Subsequently, holes are filled with particles incoming from neighboring super-cells. Remaining holes, if any, are filled using particles taken from the end of the attribute tile. In contrast, in our algorithm, particles are directly moved to their target bin-they get moved exactly once per time step.

While strict binning overcomes several of the aforementioned limitations, it is nontrivial to implement efficiently. Representing each cell by a single fixedcapacity array [START_REF] Tskhakaya | Optimization of PIC codes by improved memory management[END_REF][START_REF] Decyk | Particle-in-Cell algorithms for emerging computer architectures[END_REF] is space inefficient, and falls short outside of specific scenarios where uniform particle density can be assumed. Linked lists [START_REF] Hockney | Computer Simulation Using Particles[END_REF]Sect. 8.4.] lead to tremendous overheads both in terms of space (to represent list cells) and time (to follow indirections). Vectors, a.k.a. resizable arrays, suffer from a prohibitive 2x space overhead and involve costly resize operations.

Packed Memory Arrays (PMAs) [START_REF] Durand | A Packed Memory Array to Keep Moving Particles Sorted[END_REF] have been proposed as a specialized data structure for keeping particles sorted by their cell index. This data structure stores particles in a large array that also contains a fraction of unused cells, called gaps (a.k.a. holes). The width of the gaps may increase or decrease as particles move cells. When a gap closes, rebalancing operations involving backward or forward shifting of the particles must be performed to restore balanced gaps.

Nakashima et al. [START_REF] Nakashima | Large Scale Manycore-Aware PIC Simulation with Efficient Particle Binning[END_REF] propose a data structure that we view as a parallelismfriendly version of PMAs. To tame the frequency of rebalancing operations, the authors introduce thread-local overflow buffers. However, as the authors acknowledge [17, III.D.], these buffers come at the cost of increased complexity in the code, of additional space usage, and of slower processing of the overflow buffers, on which SIMD operations do not apply.

In Nakashima et al.'s work [START_REF] Nakashima | Large Scale Manycore-Aware PIC Simulation with Efficient Particle Binning[END_REF], most data races are eliminated thanks to the use of a 8-color scheme [START_REF] Kong | Particle-in-cell simulations with charge-conserving current deposition on graphic processing units[END_REF], which we also use. For the remaining data races, which are associated with fast-moving particles, their algorithm processes them in a separate sequential loop, which induces a major sequential bottleneck as soon as the percentage of fast-moving particles exceeds a fraction of a percent. In contrast, we are able to integrate the processing of these particles within the main parallel loop, using an atomic operation.

Furthermore, unlike prior work exploiting PMAs, our approach relies on a general-purpose bag representation, based on chunks. We only customize the bag implementation to accommodate SoA layout. Our data structure does not involve any shifting of data nor any overflow buffer. This has two main benefits. First, we save numerous memory operations. Second, the performance of our algorithm is robust to an increase in the percentage of fast-moving particles.

Future Work

In this work, we focused on multicore and SIMD parallelism. In future work, it would be great to extend our algorithm with a layer of domain decomposition, using MPI communications. We speculate that chunks could be used as buffers for emission and reception of particles reaching the cells at the frontier of a domain. These chunks could then be merged, at the end of the time step, with the locally-processed chunks. The flexibility offered by chunks might be helpful for dealing with dynamically-sized domains. Furthermore, it would be interesting to adapt our algorithm to target architectures with larger number of cores, such as GPUs or MICs. We think that the organization in chunks could help addressing the issue of load balancing, which is critical on these architectures (e.g., [START_REF] Larin | Load Balancing for Particle-in-Cell Plasma Simulation on Multicore Systems[END_REF]).

ParametersN:

  number of particles. nbCells = X • Y • Z: size of the grid. ∆t: duration of a time step. Variables particles[0..N -1]: set of particles, with position xp and velocity vp. ρ[0..X][0..Y ][0..Z]: charge density. E[0..X][0..Y ][0..Z]: electric field.

Fig. 3 .

 3 Fig. 3. Analysis of performances in the roofline model.

Fig. 4 .

 4 Fig. 4. Time evolution of electric energy in the Landau damping simulation.

Fig. 5 .

 5 Fig. 5. Time evolution of ρ(t, x, y) in the electron hole simulation, at t = 20 and t = 40.

Fig. 6 .

 6 Fig. 6. Time evolution of the y part of electric field norm in the electron hole simulation, for different values of the number of particles and of the grid size.

  8 corners per cell Foreach time step Foreach color in [0..7] // 8 coloring phases Parallel Foreach tile of that color // OpenMP parallel Foreach cell idCell in that tile Read E[x][y][z], foreach (x, y, z) among the 8 corners of cell idCell Foreach chunk in particles[idCell] Foreach particle in that chunk // SIMD vectorized Update particle velocity Foreach particle in that chunk // SIMD vectorized Update particle position Compute idCellNext, the index of the cell containing the particle Foreach particle in that chunk If the particle moves inside its tile Or it moves to the closer half of a neighbor tile Add the particle into particlesNextPrivate[idCellNext] Else Atomically add the particle into particlesNextShared[idCellNext] Add its charge into ρNext[thisCoreId][idCellNext][..] // SIMD Put a pointer to that chunk into the freelist of the current core Parallel Foreach idCell in [0..

nbCells-1] // OpenMP parallel Set particles[idCell] to particlesNextPrivate[idCell] Merge particlesNextShared[idCell] into particles[idCell] Set particlesNextPrivate[idCell] to empty, using an empty chunk Set particlesNextShared[idCell] to empty, using an empty chunk Parallel Foreach (x, y, z) in [0..X]x[0..Y]x[0..Z] // OpenMP parallel, collapsed Foreach of the 8 pairs (idCell,i) such that (x,y,z) is i-th corner of idCell Foreach idCore in [0..nbCores-1]

In practice, we allocate a dozen extra chunks per core, giving some slack and avoiding dynamic load balancing of free chunks. Note that it is very unlikely for these chunks to ever be needed, because cores free chunks at a faster rate than they fill chunks.

For a 2x2x2 tiling, at the i-th coloring phase, the algorithm processes cells whose coordinates satisfy: ((x/2) mod 2) + 2 • ((y/2) mod 2) +

• ((z/2) mod 2) = i. Using larger tiles is possible but greatly reduces the number of tiles processed in parallel.

Particle velocities in the experiment of Table1follow the sum of two Gaussian distributions, like in the bump-on-tail instability. Details may be found in[START_REF] Barsamian | Software artifacts for Euro-Par 2018 paper[END_REF].

This quantity, which we call "y part of electric field norm", is defined as half of the square root of the electric energy (E 2x + E 2 y ) minus the part of that energy corresponding to the modes in x (here, the first 20 modes).

Acknowledgments. We would like to thank the anonymous reviewers for their valuable suggestions and comments. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom Research and Training Program 2014-2018 under Grant Agreement No. 633053. Simulations were run on the EUROfusion Marconi supercomputer, in the context of the Selavlas project led by K. Kormann. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Data Availability Statement. The datasets and code generated during and/or analyzed during the current study are available in the figshare repository [START_REF] Barsamian | Software artifacts for Euro-Par 2018 paper[END_REF].


Jean-Yves Dieulot
email: jean-yves.dieulot@polytech-lille.fr

Fréderic Colas
email: frederic.colas@ensam.eu

Lamine Chalal
email: laminechalal@gmail.com

Geneviève Dauphin-Tanguy
email: genevieve.dauphin-tanguy@ec-lille.fr

Economic supervisory predictive control of a hybrid power generation plant

Keywords: Economic model predictive control, power system control, supervisory predictive control, renewable hybrid system, hardware in the loop

This work deals with the development of an economic supervisory predictive control method for the management of a hybrid renewable energy system. The hybrid cell integrates solar panels, a gas microturbine and a storage unit. Tuning the predictive controller is easy: the optimal criterion encompasses the environmental, fuel, energy delivery and storage costs. Short time predictions of the solar power are embedded in the supervisor which yields smoother battery control and better power management. Real-time experiments are driven in a Hardware-in-the-Loop framework illustrating the relevance of the proposed supervisory predictive control design.

Introduction

Renewable and alternative energy technologies are receiving worldwide attention. However, their integration in the power grid remains a challenging task. Because of changing weather (e.g. solar radiation or wind speed variations), the power of renewable sources is intermittent. One way to deal with the variable output of these energy generation systems is through the use of hybrid cells. Hybrid cells integrate renewable power sources, storage units and conventional sources in order to deliver a reference power to the grid [START_REF] Leva | Hybrid renewable energy-fuel cell system: Design and performance evaluation[END_REF][START_REF] Nehrir | A review of hybrid renewable/alternative energy systems for electric power generation: Configurations, control, and applications[END_REF][START_REF] Venkataramanan | Optimal technology selection and operation of commercial-building microgrids[END_REF]. This paper focuses on the supervision and optimal management of a hybrid power cell, which requires to consider the variety of dynamics and the versatility of the power sources. The supervisor generates the power references fed to the cell components.

In the literature, many heuristic methods have been proposed for the supervision of a hybrid cell. They use a set of rules to decide the amount of control to be applied at each time [START_REF] Boukettaya | Fuzzy logic supervisor for power control of an isolated hybrid energy production unit[END_REF][START_REF] Courtecuisse | A methodology to design a fuzzy logic based supervision of hybrid renewable energy systems[END_REF][START_REF] Hajizadeh | Intelligent power management strategy of hybrid distributed generation system[END_REF][START_REF] Hovgaard | Optimal energy consumption in refrigeration, systems -modelling and nonconvex optimisation[END_REF][START_REF] Mohamed | Real-time energy management scheme for hybrid renewable energy systems in smart grid applications[END_REF]. For example, in [START_REF] Courtecuisse | A methodology to design a fuzzy logic based supervision of hybrid renewable energy systems[END_REF], a fuzzy supervisor selects the mode of power supply and controls the load share of a distributed wind-storage power system according to the actual wind power. These methods can work model-free, but are sub-optimal. Moreover, their implementation and tuning become more complex as the number of variables increases. Global optimization algorithms have allowed to consider economic objectives. However, the intermittent sources are only seen as disturbances or uncertain negative loads and power predictions are not harnessed [START_REF] Ipsakis | Power management strategies for a stand-alone power system using renewable energy sources and hydrogen storage[END_REF][START_REF] Mohamed | System modelling and online optimal management of microgrid using mesh adaptive direct search[END_REF]. Predictive supervision strategies seem to be more suitable for hybrid cell supervision. Predictive control handles a cost function and is able to encompass predictable disturbances. Thanks to the receding horizon strategy, a predictive controller can take scheduled reference power changes or renewable power forecast into account. Model predictive control (MPC)has thus been used in several hybrid power plant control strategies (e.g. [START_REF] Guerin | Supervisory control design for systems of multiple sources of energy[END_REF][START_REF] Valenciaga | Supervisory control for a stand-alone hybrid generation system using wind and photovoltaic energy[END_REF]). There are few real-time supervisors of power systems based on economic objectives explicitly. A supervisory predictive controller with an economic criterion for the optimal management of standalone wind-solar energy generation has been proposed in [START_REF] Liu | Supervisory predictive control of standalone wind/solar energy generation systems[END_REF][START_REF] Liu | A distributed control framework for smart grid development: Energy/water system optimal operation and electric grid integration[END_REF][START_REF] Qi | Distributed supervisory predictive control of distributed wind and solar energy systems[END_REF]. However, the controller did not exactly use an Economic Predictive algorithm. Indeed, artificial weights were introduced in the quadratic criterion in order to smooth the battery State Of Charge (SOC). In [START_REF] Hovgaard | The potential of economic mpc for power management[END_REF], a simple linear economic criterion was introduced for power flow management. The non-convex Economic MPC control of a refrigeration unit, which included energetic considerations, was proposed in [START_REF] Hovgaard | Optimal energy consumption in refrigeration, systems -modelling and nonconvex optimisation[END_REF]. This paper fills in a gap and proposes a true economic predictive supervisor for a hybrid power cell. Contrary to existing predictive supervisors, the economic criterion is only based on actual costs and is a nonstandard function -i.e. not as for a standard Linear Quadratic regulator -of the control variables. Moreover, the supervisor embeds explicitly the renewable power predictions. The objective of the present work is to develop an Economic Supervisory Predictive controller (ESP) for the optimal operation and management of a hybrid cell. This cell comprises a microturbine, a storage unit and solar panels. The ESP will compute the reference powers according to financial considerations, encompassing renewable energy tariffs as well as storage, fuel and environmental costs. At first, a closed-loop model of each component of the hybrid cell will be presented. Next, the design of the objective function will be addressed and the predictive control problem will be solved. Eventually, real-time experiments in a Hardware In the Loop (HIL) framework, i.e. mixing real and emulated components, will show the relevance of the economic supervisory predictive approach.

HYBRID CELL ARCHITECTURE AND COMPONENT CON-TROL

The different components of the hybrid cell (an emulated microturbine, real solar panels and a real storage unit) are connected to the grid via electronic converters (Fig. 1). Since the microturbine and solar panels are fixed, the sizing of the storage unit and the power planning are subsequently designed using the software HOMER. The optimization is done on the basis of a one year range. The time-scale of the hybrid cell supervision is much shorter, as the prediction horizon value is only 10 s. In the framework of the design and the validation of the predictive supervisory controller, three kinds of models need to be considered:

• a general (closed-loop) model for simulation purposes,

• a simplified model, if possible linear, in order to design the supervisor,

• real-time devices or emulated models which are used for validation.

The predictive supervisor will, in this case, deliver the power references to the storage unit and the microturbine. The dynamical models, which are described thereafter, will be embedded into the optimization problem.

Solar power prediction

The photovoltaic system consists of 108 modules BP solar 3160; each delivers a power of 160W. These modules are connected to a 3-phase grid via an inverter. As will be shown later, the predictive supervisor can take advantage of solar power short-term predictions. Numerous solar power forecasting methods have been proposed so far, such as physical and statistical approaches, parametric models and trend curves, etc. [START_REF] Paolia | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF]. Out of several trend functions (linear, quadratic, exponential), the linear trend gives the best results. Only the comparison between the estimated and the measured values is shown ( Fig. 2). The solar power is measured during a cloudy day. The prediction horizon H p is equal to 10 s. This predictor is simple and accurate and thus fits for real-time implementation. As will be seen later, a solar power profile can be reused (stored and then emulated) in the HIL framework to compare several algorithms.

Storage unit

The storage unit combined with the microturbine is used as a backup for the solar array. The model describing the relationship between the voltage, current and the SOC of a battery can be found in [START_REF] Chungpaibulpatana | A simulation model for predicting the performance of a solar photovoltaic system with alternating current loads[END_REF]. A simplified model is used to represent the battery storage:

P bat = 1 τ p s + 1 Pbatref , (1)
Pbatref = max(min(P batref , P max ), P min ),

where s is the Laplace operator, P batref , P bat are the storage unit reference and real output powers; Pbatref is the saturated reference power and P max , P min are respectively the maximum power of charge and the maximum power of discharge.

It is assumed that the storage unit SOC is a function of the integral of the power (in the working range). In the case of SOC saturation, if t 0 P bat dt > E max or t 0 P bat dt < E min , then P bat = 0. E max and E min are respectively the maximum and minimum stored energy. τ p = 5s is the time constant of the battery.

When there is no saturation, one has a first-order equation:

P bat = 1 τ p s + 1 P batref . ( 3
)
This model is implemented on a dedicated test bench which comprises a set of supercapacitors. Fig. 3 shows the principle of this storage unit, for which the power reference can be generated through a power amplifier. The supercapacitors are connected to the grid via two converters, an interlaced chopper (reduces current ripple) and an inverter. The supercapacitors consist of 6 Maxwell modules in series, the characteristics of each being 48V, 160F. In this case, the models used for design and validation are identical. The charge and discharge of the storage unit is fast, which is a limitation with respect to a lead-acid battery but will enlighten the relevance of the SOC management.

Microturbine generator

Microturbine generators are known to deliver clean energy from a wide range of fuels with a low level of gas emissions. As traditional gas turbines, they have a gas compressor, a gas turbine and a recuperator, along with a high speed permanent magnet synchronous generator, power electronic devices and their grid interface. In this paper, a simplified model is taken into account without recuperator (considered as a simple heat exchanger). Also, complexities with start-up and shut-down operations are discarded.

Generally, the control structure of a MTG includes three important control blocks monitoring respectively speed, acceleration and temperature [START_REF] Pai | A design and operation of power converter for microturbine powered distributed generator with capacity expansion capability[END_REF][START_REF] Zhu | Development of models for analyzing the loadfollowing performance of microturbines and fuel cells[END_REF]. The main control unit consists of the speed controller. The acceleration controller is mainly used to limit the rate of rotor acceleration at start-up. The temperature controller restricts the temperature at a predetermined firing temperature. The simplified model only retains the speed loop, as shown in Fig. 4. One has:

C m = 1.3(W f -0.23) + 0.5(1 -N ), (4)
N = 1 T l s + 1 C m , (5)
W f = 1 T f s + 1 a cs + b (0.77N.V CE + 0.23), (6)
The speed controller is:

V CE = 25 0.05s + 1 (N ref -N ), (7)
where C m , N, W f are the per unit rotor torque, rotor speed, fuel flowrate, V CE is the per unit fuel command per unit speed, the turbine power P mtg is the product of speed and torque. T l , T f , a, b, c are respectively the turbine rotor time constant, fuel system time constant, and fuel system transfer function parameters. This model is implemented into the real-time simulator for validation. In the case of "stiff" turbines, it is possible to simplify the closed-loop model of the microturbine, which finally reduces to a fourth-order transfer function:

P mtg = 1 (T 1 s + 1)(T 2 s + 1)(T 3 s 2 + T 4 s + 1) P mtgref , (8)
where T 1 = 0.05, T 2 = 0.15, T 3 = 0.22 and T 4 = 0.27. This coarse linear model shows a good agreement with the full model in the region where the microturbine operates when embedded in the hybrid cell. This model is used in the Model Predictive controller design procedure. One can see in Fig. 5 Figure 5: Actual and reference powers that the electrical power produced by the microturbine follows the required reference; simulations are consistent with that found in the example by [START_REF] Martinez | Tools for analysis and design of distributed resourcespart iv: Future trends[END_REF].

Finally, the simplified dynamical model of the cell along with the storage unit and the microturbine can be put under the discrete-time state-space form:

x k+1 = Ax k + Bu k , y k = Cx k (9)
where the subscript k indicates the k th sample , x is the state space, the vector of controllable powers, y = (P bat P mtg ) T is the output, and the vector of power references u = (P batref P mtgref ) T is the input vector that will be computed by the supervisor.

DESIGN OF AN ECONOMIC PREDICTIVE SUPERVISOR

Design of an economic criterion

Model Predictive Control (MPC) is a finite horizon optimization algorithm applied to a discrete-time model of a dynamical system. A controller which minimizes a cost function depending on the state and input of the process is computed over a time window [t, t + H p ]. The possible state trajectories that can be obtained from the current state are explored using a prediction based on the dynamical model. Via Euler-Lagrange equations, a sequence of optimal controls can be found which minimizes the optimization cost function over the prediction horizon. Typically, in our case, the criterion is a quadratic function of the future controls, the use of Sequential Quadratic Programming (SQP) optimization packages is used which provides, at every instant, the sequence of relevant control values. The first step of the control strategy is implemented, the calculations are carried out from the new value of the discrete-time state. Hence, the MPC is called a receding horizon strategy as the horizon is shifted forward. Hence, the main objective of the ESP controller is to manage the power flow of the hybrid cell by minimizing the actual operational cost over a short-term receding horizon. The costs considered are the price for gas consumption C fuel , the tax imposed by the government on carbon dioxide emissions C emissions , the battery cycling cost C cycling , and a cost that accounts for tariff policies. The controller determines the power reference of the microturbine P mtgref and the power reference of the storage system P batref at each time-step over a finite horizon. One has to solve an optimization problem that minimizes the cost function subjected to model and operational constraints. The microturbine operates in closed-loop with the local speed controller described previously.

The prediction horizon H p is chosen by a tradeoff between solar prediction accuracy (which decreases with the prediction horizon length) and the other components' time constants. In this work, the prediction and control (that is, optimization) horizons are the same. Fig. 6 shows the simplified control scheme. The total power delivered by the cell P t is :

P t = P pv + P bat + P mtg . (10)
The discretized performance index at control step k over a prediction horizon H p is as follows:

J = Hp k=0 α(P d,k -P pv,k -P bat,k -P mtg,k ) 2 + Hp k=0 C fuel,k + Hp k=0 C emissions,k + Hp k=0 C cycling,k (11)
where the subscript k corresponds to the value of a variable at sample k. The financial cost J can be expressed in e. The weight α ((e/W 2 ) is not artificial but corresponds to a tariff policy defined by the grid manager which penalizes the gap between the required power P d,k and the delivered power P t,k = P pv,k + P bat,k + P mtg,k . P bat,k , P mtg,k correspond respectively to the storage unit and microturbine powers to be determined by the algorithm. P pv,k corresponds to the solar power predictions over the receding horizon. The reference power delivered by the storage unit and the microturbine will be equal to P d,k -P pv,k . Hence, the accuracy of the predictions is quite important. More complications can be added in the same principle of ESP, e.g. penalty criterion from the grid manager. The individual contributions are given thereafter.

A-Fuel cost model Since little is known about the operating cost of a real microturbine, it is chosen to use the cost function of a diesel generator, C fuel . This cost is expressed as a non-linear function of the produced power P mtg [START_REF] Dufo-Lopez | Design and control strategies of pvdiesel systems using genetic algorithms[END_REF] and is shown in the following equation:

C fuel = C F [BP Nmtg + AP mtg ], (12)
where :

• C F is the fuel price (e/l),

• A = 0.246 l/kW and B = 0.08415 l/kW are the fuel curve coefficients,

• P Nmtg is the rated power (kW),

• P mtg is the generator output power(kW),

• C fuel is the fuel consumption cost (e).

The cost is thus an affine function of the microturbine power.

B-Pollutant emissions (CO 2 )

CO 2 emission E is directly related to the fuel consumption, which can be approximated by a quadratic function of the generator power P mtg as [START_REF] Johansen | A qualitative approach to economic-environmental dispatch[END_REF] :

E(P mtg ) = A + BP mtg + CP 2 mtg (13)
where A, B, C are scalar coefficients. As an example, France issued a carbon tax (17 e/ton) with a new levy on oil, gas and coal consumption:

C emisions = 0.017 * E(P mtg ). (14)

C-Battery cycling cost

Battery cycle lifespans vary according to their type and monitoring policy. The value of a mean cycle (full charge and full discharge as shown in Fig. 7) is given by the total cost of the storage unit divided by the lifespan (in cycles): cycle cost = cost of the battery number of expected cycles [START_REF] Qi | Distributed supervisory predictive control of distributed wind and solar energy systems[END_REF] Since the battery usually does not finish a complete cycle, the cost of a portion of cycle is given by:

portion of a cycle cost = 1 2 E B -E A E A (cycle cost) (16)
The factor 1/2 accounts for the charge and discharge of the battery during a full cyle (see Fig. 7 ). Assuming that the SOC varies linearly, the cycling cost can be defined as: where λ can be calculated from equation [START_REF] Hovgaard | The potential of economic mpc for power management[END_REF].

C cycling = λ dSOC dt , (17)
Traditional predictive controllers often use explicit quadratic criteria which involve the square of the control input and the square of the output error. The criterion presented in this paper has also linear and quadratic terms, but these are obtained in a quite different way, only on an economic basis.

Working out the predictive controller

As is said before, the predictive control strategy computes the storage unit and the microturbine power references by minimizing the economic criterion J. Additional operational constraints which limit the reference power change rate [START_REF] Liu | Supervisory predictive control of standalone wind/solar energy generation systems[END_REF] need to be considered. The microturbine must work at a load level more than 50% of its nominal power, let otherwise the efficiency will decrease drastically:

15kW ≤ P mtgref ≤ 30kW. ( 18
)
Inner-loop constraints on the fuel flowrate F mtg read:

0 ≤ F mtg ≤ F max . ( 19
)
The battery operational constraints include power and capacity limitations: improper depletion or battery overcharge should be avoided. The battery SOC is kept within reasonable limits to increase life duration of the battery:

-5kW ≤ P batref ≤ 5kW. ( 20)

0.5 ≤ SOC ≤ 0.8. ( 21
)
and

dSOC dt ≤ 0.1. ( 22
)
Instead of imposing a maximum rate of variation for the reference or real powers, it has been chosen to restrain the variation of the SOC. It is not mandatory to restrain formally the variation of the microturbine power which is already limited by a comparatively slow dynamics. Eventually, it is not possible to enforce the rate of variation of the total power which embeds the possibly fast variations of the solar power. The overall predictive control problem can be put under the following form, at each sample time: minimize

P batref,k ,P mtgref,k J subject to ( 18)( 20)( 21)( 22)( 19), ∀k ≤ H p . x k+1 = Ax k + Bu k , y k = Cx k (23)
As was shown in the previous subsection, the criterion depends explicitly, in a Linear Quadratic form, on the storage unit and microturbine powers P bat,k and P mtg,k . The computation of a solution to the QP problem is handled using a standard optimization routine. A sequence of reference trajectories is found, only the first value is used by the ESP, and the procedure is repeated at each sample time.

HARDWARE IN THE LOOP VALIDATION

Power Hardware In the Loop Structure

The Power HIL principle is applied on the RT-LAB simulator. This simulator can run real time models on a multi-CPU computer and emulate some physical components of the hybrid renewable power plant (for more details about implementation, see [START_REF] Martinez | Tools for analysis and design of distributed resourcespart iv: Future trends[END_REF]). The architecture used in this paper is illustrated in Fig. 8. The power amplifier interconnects the simulated components to the real ones. Hence, it is possible to create a real power from a (simulated) component model through the real-time simulator. The same is done for the microturbine represented previously in Fig. 4. This simulated source is blended with the real power measurements of the solar panels and of the real storage unit in order to deliver the required total cell power. Standard solar power profiles are used to make a fair comparison of algorithms. Power measurements from solar panels are collected for the first time and then emulated with the real time simulator, ensuring repeatable external conditions. In our case, the ESP is designed with a receding horizon H p = 10s, references are sent every T = 1s, and the sampling time is fixed to 0.01s.

Real-time simulations

At first, some simulations are undertaken to verify the feasability of the approach. Table 1 is showing the good match between simulations and realtime results for a standard profile. The Mean Relative Absolute Prediction Error (MRAPE) for each of the simulated and real components' powers never exceeds 10 %. The MRAPE for the total power is below 2 % . During the first real-time experiment, the reference power is kept to 30 kW and the solar power is slightly variable. In order to enlighten the relevance of integrating the solar predictions into the supervisor, results are compared with the same controller working without predictions. In the latter case, the solar power is considered as a measured disturbance as in [START_REF] Liu | Supervisory predictive control of standalone wind/solar energy generation systems[END_REF]. Results are also compared with that given by a rule-based controller for which the decision tree is given in Fig. 9 (only for off-line simulations). As expected, the predictive controllers outperform a basic rule-based controller. Results are shown in Table 2. The Mean Relative Absolute Error between the total power and the power reference when using the controller proposed in this paper is M RAE = 1.86%. The microturbine and battery powers are maintained within their respective limits. The use of predictions does not improve the total cost but allows a better management of the battery power. Indeed, the battery power profile is more smooth when predictions are used which is confirmed by the evolution of the SOC shown in Fig. 10. Fig. 11 and Fig. 12 are showing the experimental results respectively for a standard solar profile when the solar power is predicted and not predicted. When the solar power is not predicted, the storage unit compensates the fluctuations of the solar power directly. However, when solar power predictions are used, the power tracking is worse because of the discrepancy between real and predicted solar powers. Since the total cost is almost the same, one will prefer to use solar power predictions in order to obtain better battery management, and leave them when better power tracking is desired.

Real-time experiments for an ESP with predictions are now proposed for a steady and an unsteady solar power profiles in order to evaluate the behavior of the algorithm. The following figures are showing the corresponding powers and SOC. When the solar power is steady, the variations of the microturbine and battery powers are small, , according to the results in Fig. 13. The drift in the SOC of the storage unit can reach 5% as seen in Fig. The tracking of the total power reference is nearly perfect, and, with respect to a one-step-ahead algorithm, the results will be little improved. On the contrary, the predictive approach is more interesting when the solar power 16 that the storage unit is quickly used and then hands over to the microturbine, which dynamics is much slower. The microturbine will in turn allow the battery to recharge while keeping the total power close to its reference value (see 15). This behavior is not really intuitive; tracking is not very good when the solar power varies quickly due to the components' inertia and solar power prediction errors . Moreover, it could be interesting to see if better tracking and smoothing of the power profiles could be achieved. Instead of enforcing hard constraints on the power rates of variation, an alternative would consist of detecting abrupt changes, then decrease the value of the receding horizon. This will give a more reactive control during a short time. On the contrary, increase the receding horizon for steady operating conditions.

Conclusion

The economic supervisory predictive control of a hybrid renewable energy system is designed according to technical and financial considerations. The supervisor sends the reference power to every controllable component of the cell. The design of a nonstandard quadratic and economic criterion and the use of dynamical models allow rapid and effective tuning of the supervisor. The predictive strategy allows to take short-time renewable power predic- Hence, this supervision strategy can be implemented in real-time and with true components. The experimental results are close to that expected from off-line simulations and show the relevance of the theoretical developments. Several solar power profiles have been tested, showing that results are quite good when the solar power varies in reasonable ranges. When variations are much faster, the tracking is not so good, because of the possible discrepancies between the real and predicted solar powers. It can be expected that an adaptive algorithm where the prediction horizon is a function of the solar power prediction accuracy, that is, the horizon is shrunk when fast solar power variations occur, would improve the existing results. Such an adaptive approach is the topic of future research. Another perspective is to generalize the methodology to large-scale units, by incorporating new elements.

Figure 1 :

 1 Figure 1: Scheme of the hybrid renewable power plant

Figure 2 :

 2 Figure 2: Measured and predicted photovoltaic power

Figure 3 :

 3 Figure 3: Storage unit test bench

Figure 4 :

 4 Figure 4: Main block of microturbine model adapted from[START_REF] Pai | A design and operation of power converter for microturbine powered distributed generator with capacity expansion capability[END_REF]

Figure 6 :

 6 Figure 6: Simplified control scheme

Figure 7 :

 7 Figure 7: Full discharge and full charge of a battery

Figure 8 :

 8 Figure 8: Power Hardware In the Loop architecture

  Fig. 14.

Figure 9 :

 9 Figure 9: Flowchart of a rule-based controller

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: SOC with and without taking into account the solar predictions

Figure 13 :Figure 14 :Figure 15 :Figure 16 :

 13141516 Figure 13: Power of each element for a steady profile

Table 1 :

 1 Mean relative absolute error "%"

	Total power Turbine power Battery power
	1.7	6.4	8.3

Table 2 :

 2 Cost function components

		hierarchical ESP without ESP with
		algorithm	predictions predictions
	Emission	358	385.58	371.82
	CO2 (Kg)
	Emission	6.1	6.55	6.32
	cost (e)
	Fuel cost	55.33	30.14	31.45
	(e)
	Battery	1.78	5.83	4.76
	cycling
	cost (e)
	Total cost	63	42.52	42.53
	(e)

varies more significantly. One can see from Fig.