Simon Bernard
email: simon.bernard@univ-lehavre.fr

Feriel Chikh-Bled

Hasna Kourchi

Farid Chati

Fernand Léon

Broadband negative reflection of underwater acoustic waves from a simple metagrating : modeling and experiment

Metagratings are periodic arrays of subwavelength scatterers, or atoms, engineered to refract or reflect waves towards anomalous directions with unit efficiency. Here we design and build a metagrating to control the reflection direction of water-borne ultrasound waves impinging on rigid or free surface. The grating and its atoms are designed to cancel the specular reflection and to redirect acoustic power towards a negative reflection direction, through the first negative order Floquet mode. Despite a simple design, based on c-shaped brass particles acting as Helmoltz resonators, the grating is efficient (> 90%) over a relatively broad range of frequencies (74-103 kHz) for a broad range of incident angles (14 to 54 • ). This good performance is obtained by tuning the distance between the atoms and the reflective surface. A multiple scattering analytical model is presented to explain the phenomenon, and a finite element model is developed to further investigate the performance of the proposed design. Predictions from the model are confirmed experimentally in a water tank. The simplicity, reconfigurability and scalability of the design, as well as its high efficiency, broadband behavior and robustness to incident angle are all features that make the grating potentially useful for various applications in underwater acoustic such as telemetry, communication, or noise mitigation.

I. INTRODUCTION

Bulk 3D metamaterials and thin 2D metasurfaces are designed to offer extraordinary wave propagation properties, usually from an assembly of multiple subwavelength elements. Many metasurfaces have been developed in the past years, both for electromagnetic waves [START_REF] Holloway | An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials[END_REF][START_REF] Pfeiffer | Metamaterial Huygens' surfaces: tailoring wave fronts with reflectionless sheets[END_REF][START_REF] Yu | Flat optics with designer metasurfaces[END_REF][START_REF] Zhao | Recent advances on optical metasurfaces[END_REF][START_REF] Glybovski | Metasurfaces: From microwaves to visible[END_REF] and acoustic waves [START_REF] Ma | Acoustic metamaterials: From local resonances to broad horizons[END_REF][START_REF] Assouar | Acoustic metasurfaces[END_REF][START_REF] Liang | Wavefront manipulation by acoustic metasurfaces: from physics and applications[END_REF][START_REF] Wu | Perspective: Acoustic metamaterials in transition[END_REF]. One application of metasurfaces is to redirect waves toward anomalous directions not classically allowed by the reflection and refraction laws [START_REF] Li | Reflected wavefront manipulation based on ultrathin planar acoustic metasurfaces[END_REF][START_REF] Li | Experimental realization of full control of reflected waves with subwavelength acoustic metasurfaces[END_REF][START_REF] Mei | Controllable transmission and total reflection through an impedance-matched acoustic metasurface[END_REF][START_REF] Xie | Wavefront modulation and subwavelength diffractive acoustics with an acoustic metasurface[END_REF]. This can be achieved using a phase-gradient metasurface : a local phase shift is applied by the surface to modify the wave's momentum and reroute it toward the desired direction, according to the so-called generalized Snell-Descartes's laws [START_REF] Yu | Light propagation with phase discontinuities: generalized laws of reflection and refraction[END_REF]. The phase gradient can be discretized by engineering multiple small elements, or atoms, to offer the desired local phase shift. In acoustics, space-coiling structures or Helmholtz-like resonators are often employed [START_REF] Assouar | Acoustic metasurfaces[END_REF], but require tedious design and precise manufacturing. In air acoustics, 3D printing can provide complex structures quickly and at low cost, but most 3D printing materials are inappropriate for underwater applications because of a low contrast in acoustic impedance. Furthermore, such surfaces usually work optimally in a narrow frequency band, as complex meta-atoms will offer the desired phase response only at their design wavelength. Recent works have highlighted another limitation : even with a perfect realization of the phase gradient, the efficiency is intrinsically limited for large steering angles [START_REF] Estakhri | Wave-front transformation with gradient metasurfaces[END_REF][START_REF] Díaz-Rubio | Acoustic metasurfaces for scattering-free anomalous reflection and refraction[END_REF]. This can be overcome by designing locally active or non-local metasurfaces, relying on near field coupling between the atoms [START_REF] Díaz-Rubio | Acoustic metasurfaces for scattering-free anomalous reflection and refraction[END_REF][START_REF] Epstein | Synthesis of passive lossless metasurfaces using auxiliary fields for reflectionless beam splitting and perfect reflection[END_REF]. However, this comes at the cost of further increased design and fabrication challenges.

The concept of metagrating relies on a different physical principle and on a simpler design based on a periodic cell containing down to one single atom [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Epstein | Unveiling the properties of metagratings via a detailed analytical model for synthesis and analysis[END_REF][START_REF] Torrent | Acoustic anomalous reflectors based on diffraction grating engineering[END_REF][START_REF] Hou | Highly efficient acoustic metagrating with strongly coupled surface grooves[END_REF][START_REF] Fu | Multifunctional reflection in acoustic metagratings with simplified design[END_REF][START_REF] Ni | High-efficiency anomalous splitter by acoustic metagrating[END_REF][START_REF] Fan | Metagratings for waterborne sound: Various functionalities enabled by an efficient inverse-design approach[END_REF][START_REF] Chiang | Reconfigurable acoustic metagrating for high-efficiency anomalous reflection[END_REF]. The possible outgoing wave directions are governed by the periodicity (Fig. 1). By engineering the atoms, the unwanted diffraction channels can be suppressed and all the impinging energy rerouted to a desired direction with unitary efficiency [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Torrent | Acoustic anomalous reflectors based on diffraction grating engineering[END_REF]. Such gratings have been recently demonstrated experimentally for acoustics in air [START_REF] Chiang | Reconfigurable acoustic metagrating for high-efficiency anomalous reflection[END_REF] and numerically for underwater acoustics [START_REF] Fan | Metagratings for waterborne sound: Various functionalities enabled by an efficient inverse-design approach[END_REF]. Another benefit over phasegradient metasurfaces is the existence of configurations effective over a large bandwidth and a large range of incident angles, as recently demonstrated in optics [START_REF] Deng | A facile grating approach towards broadband, wide-angle and high-efficiency holographic metasurfaces[END_REF][START_REF] Estakhri | Visible light, wide-angle graded metasurface for back reflection[END_REF][START_REF] Deng | Facile metagrating holograms with broadband and extreme angle tolerance[END_REF] and in this study for acoustics.

Resonant atoms can provide the desired scattering properties [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF]. Here we use c-shaped atoms, which behave as Helmholtz resonators with a low frequency resonance mainly governed by the size of the cavity and neck opening. Such atoms have recently attracted interest in acoustic wave manipulation, as they are easy to fabricate and offer interesting scattering properties [START_REF] Chiang | Reconfigurable acoustic metagrating for high-efficiency anomalous reflection[END_REF][START_REF] Melnikov | Acoustic metamaterial capsule for reduction of stage machinery noise[END_REF]. In this study, we design and build a metagrating for anomalous reflection of underwater ultrasound waves around 90 kHz, from small (radius ≈ λ /6) brass atoms, with a focus on negative reflection : the goal is to redirect the incoming energy toward the direction of the incident wave. We demonstrate how, by adequately choosing the distance between the grating and the reflective surface, it is possible to obtain near unit efficiency over a large frequency band and range of incident angles.

The broadband and broad-angle effects are linked to the multiple scattering of the wave by the surface and the atoms, and an analytical model is provided, as well as Finite Element Modeling and an experimental demonstration.

z=0 ¡ r (deg) m = -1 m = 0 p m = 0 m = -1 m = 1 p = ¢ ¡ i (deg)
retro-re ection :

¡ i = -¡ r = 30° ¡ r = asin(m¢/p+sin ¡ i ) (a) Metagrating (b)
Re ective boundary (rigid or free)

¡ i h ¡ r x z x=0 FIG. 1. (a)
To achieve anomalous reflection, metagratings rely on periodicity, which allows only a finite number of reflected modes to exist, given by the grating equation, illustrated in (b) for p = λ . By tailoring the scattering properties of the atoms, the m = 0 mode and the specular reflection from the surface behind the grating can be tuned to cancel each other so that all the energy is rerouted toward the desired anomalous direction.

II. ANALYTICAL MODEL AND PARAMETER SELECTION

Consider a periodic grating of small atoms located at a height h from a perfectly reflecting surface. The surface can be either an acoustic free (acoustic pressure = 0) or a rigid (normal velocity = 0) boundary condition. When the grating is illuminated by a plane acoustic wave with incidence angle θ i , only a finite number of reflected plane waves, or Floquet modes, can exist.

Their reflection angles θ r are related to the periodicity p and wavelength λ by the grating equation [START_REF] Popov | Introduction to diffraction gratings: Summary of applications[END_REF][START_REF] Maystre | Analytic properties of diffraction gratings[END_REF]:

sin(θ r ) -sin(θ i ) = m λ p , (1)
where m is an integer. Solutions of Eq. ( 1) are plotted on Fig. 1(b) for p = λ . In this case, three reflection channels exist : m = 0, for which θ r = θ i for all θ i , and the m = -1 and m = 1 modes, existing respectively for θ i > 0 and θ i < 0. With these parameters, θ i = 30 • corresponds to the case of retro-reflection, where the m = -1 mode is reflected exactly toward the direction of the incoming wave. When two or more channels co-exist at a given incidence angle, their relative weight depends on the scattering properties of the atoms [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF]. The goal is to tune the scattering such that the weight of one selected mode is made as close as possible to 1, while other unwanted modes are close to 0. In this study we focus on negative reflection : in the range of parameters for which the m = 0 and m = -1 modes exists, we design the grating such that the m = 0 mode cancels the specular reflection from the surface behind the grating, so that all energy is reflected by the m = -1 channel, in a negative direction as shown in Fig. 1(a).

Under a monopolar scattering approximation, valid for small atoms, the m = 0 mode exactly cancels the specular reflection if, for the rigid boundary condition, (see the Appendix)

1 α = 2i p ωc cos θ i cos 2 (kh cos θ i ), (2)
and for the free boundary condition

1 α = -2i p ωc cos θ i sin 2 (kh cos θ i ), (3)
where α is the effective monopolar polarizability of an atom and k = 2π/λ = 2π f /c. We see in both cases that α must be purely imaginary, which happens at the metagrating resonance [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Tretyakov | Analytical modeling in applied electromagnetics[END_REF].

Note that α is the effective polarizability of an atom, defined in the presence of all other atoms and reflective surface. Futhermore, in order to perfectly match the power of the incident wave to power radiated through the m = -1 channel we need to satisfy [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF] (see also the Appendix) :

1 cos θ i cos 2 (kh cos θ i ) = 1 cos θ r cos 2 (kh cos θ r ), (4)
for the rigid boundary condition, and

1 cos θ i sin 2 (kh cos θ i ) = 1 cos θ r sin 2 (kh cos θ r ), (5)
for the free boundary condition. We see, as noted by Ra'di et al. [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF], that retro-reflection (θ r = -θ i ) is a trivial solution of these equations, meaning that it can be obtained whatever the distance h. However, as we will show now, the solution is more robust to a change in wavelength or angle for some values of h, providing a broadband and broad-angle behavior.

We define ε as the ratio of left hand to right hand term in eq. ( 4) or [START_REF] Glybovski | Metasurfaces: From microwaves to visible[END_REF]. By considering wave scattering by the grating and its virtual conterpart on the other side of the reflective surface, we can arrive at (see the appendix) ε is the ratio of of left hand to right hand term in eq. ( 4) or [START_REF] Glybovski | Metasurfaces: From microwaves to visible[END_REF]. A value of 1 (white color) satisfies the conditions for anomalous reflection with unitary efficiency. The vertical white line at 87 kHz is the trivial solution θ r = -θ i . For h = 8 mm for the rigid surface, ε is close to one over a broad range of frequencies, particularly in the range [80-105] kHz. For the free surface, two interesting solutions exist : h = 13 mm and h close to 0 (the smallest possible value in practice being equal to the atom's radius).

ε = cos θ r cos θ i cos 2 (kh cos θ i ) cos 2 (kh cos θ r ) , (6)
for the rigid boundary condition, and

ε = cos θ r cos θ i sin 2 (kh cos θ i ) sin 2 (kh cos θ r ) , (7)
for the acoustic free boundary. We plot ε as a function of h and incident wave frequency f in Fig. 2 for both boundary conditions, for θ i = 30 • , c = 1470 m/s, and p = λ = 17 mm (and using the grating equation Eq. ( 1) to obtain θ r as a function of f ). In addition to the trivial solution (θ r = -θ i ) corresponding to the vertical white lines at about 87 kHz visible in Fig. 2(a)-(b), other solution branches exist. In the rigid case, for h = 8 mm we observe a large valley of near unity ε, due to a solution branch crossing the θ r = -θ i solution. In this region, the metagrating is able to efficiently couple the incoming power to the m = -1 channel over a broad range of frequencies, through complex wave interactions between the grating and the reflective surface. In the free boundary case, two interesting values of h exist : h very small, with the minimal values in practice being equal to the radius of the atoms (2.5 mm here, see below), or h = 13 mm.

The dispersion relation between the wavelength and reflective angle (Eq. 1) implies that, if the specular reflection is canceled over a broad frequency range, the energy is reflected toward a varying angle. Consequently, a broadband wave cannot be redirected to a given direction. This is a limitation of the metagrating. On the other hand, we can also observe that reflective and incident angle have a symmetric role in Eqs. ( 4)-( 5), implying that the broadband behavior is associated to robustness with respect to the incident angle.

III. FINITE ELEMENT MODELING OF THE PERIODIC GRATING

A Finite Element Model (FEM) of the periodic metagrating was built using Comsol Multiphysics. The acoustic structure interaction module was used to take into account the elastic nature of the atoms. Geometry and material properties were set according to the actual grating (see Similar results are displayed on Fig. 4 for the free boundary case, for height h = 13 mm and h = 2.5 mm respectively on panels (a-b) and (c-d). For h = 13 mm, the observed effect is qualitatively similar to the effect observed for h = 8mm for the rigid boundary. This is not surprising, as in both case the broadband effect is due to a crossing of two solution branches of the equation ε = 0 in the frequency -distance h space (Fig. 2). For h = 2.5 mm the effect is qualitatively different : the range of incidence angle over which the efficiency is > 0.9 is even larger, ranging from 10 to 78 • .

The grating is therefore highly efficient from near normal to very grazing incidence. Frequency band of high efficiency is somewhat narrower in this case, but the decrease in efficiency when moving away from the design frequency is less sharp than in the previous cases and the grating still has > 0.5 efficiency on a reasonably broad band (67 to 107 kHz at 30 • incidence).

IV. EXPERIMENTAL DEMONSTRATION WITH A RIGID SURFACE

Experimentally, the atoms making up the metagrating were made from 300 mm long brass hollow cylinders, with external radius r ext = 2.50 mm and internal radius r in = 2.05 mm. They were opened on one side using a milling cutter to form c-shaped atoms (see Fig. 5). The opening was δ = 2 mm large, forming a Helmholtz resonator with a broad resonance at about 90 kHz. Brass was selected because 1) it could be found easily in the desired shape at low cost, 2) it is resistant to corrosion, even when immersed in water for several consecutive days, 3) it can be easily milled (as compared to stainless steel), and 4) it is stiff (E = 100 GPa) and heavy (ρ = 8500kg/m 3 ), providing a large acoustic impedance contrast with water. A pairs of PVC blocks were drilled to the desired periodicity p = 17 mm (≈ λ at 87 kHz) to assemble the grating from 17 atoms.

the structure was then attached in front of a 10 mm thick stainless steel plate (L = 300 by 300 mm) acting as an almost perfect reflector for ultrasound waves in the frequency range of interest.

The distance between the center of the tubes and the plate was set to h = 8 mm, according to the analysis above.

Experiments were conducted in a water tank (3 m diameter, 1.5 m height). Two identical ultra- is presented on Fig. 5. A pulser-receiver (5052R, Panametrics) was used to generate a broadband input signal. The received signal was recorded on a digital oscilloscope (HDO6034, Teledyne Lecroy). The receiving transducer was rotated over 180 • in front of the plate, with 2 • steps. For each receiving position, 100 recordings were averaged to improve signal-to-noise ratio. Time windowing was applied to isolate the reflection on the grating from reflections on the tank boundary and water surface. The spectrum was obtained from Fast Fourier Transform and frequency band correction was applied using the spectrum of the specular reflection on the steel plate alone. The experiment was repeated for three incidence angles (15

• , 30 • , 45 • ).
The experimentally recorded reflected pressure is displayed as a function of the receiving angle and frequency on Fig. 6 We can also observe narrow gaps in the reflected mode spectra, at frequencies corresponding to the appearance of an additional mode, which are the acoustic analog of the Wood's anomaly [START_REF] Liu | Investigation of the origin of acoustic wood anomaly[END_REF][START_REF] Maradudin | Rayleigh and wood anomalies in the diffraction of acoustic waves from the periodically corrugated surface of an elastic medium[END_REF], whose discussion is outside the scope of the present work. with the metagrating for incidence angles 15 • , 30 • , and 45 • , respectively. The negative reflection effect is clearly visible, as in all three cases with the grating, most of the wave pressure is reflected backward, towards the quadrant containing the incident wave vector. The pressure amplitude increases compared to specular reflection pressure for large negative reflection angles (panel (b)), and decreases for small negative reflection angles (panel (d)). This is expected, as a consequence of energy balance: the power carried away by the anomalous reflection mode is proportional to the product of acoustic pressure and cosine of the reflection angle (see eqs. A.16 and A.17 in the Appendix). For 30 • incidence angle (panel (c)), the amplitude should be expected to be equal to the specular reflection pressure (panel (a)), but is slightly lower. This is an artefact due to our experimental system : at -30 • reflection angle (perfect retroreflection) the receiver is located in front of the acoustic source, actually reducing the amount of incident energy.

θ i =15°( c) θ i =30°(d) θ i =45°( e) θ i =15°( f) θ i =30°( g) θ i =45°θ r (deg)

V. DISCUSSION AND CONCLUSION

In this work, we designed a metagrating for negative reflection of underwater acoustic waves.

This metagrating can operate in front of either a rigid or an acoustic free surface. A detailed analysis of multiple wave scattering between the grating and the surface was presented. This analysis showed that the height h of the grating above the surface is a critical parameter regarding performances. Finite element modeling confirmed that, once an appropriate value for h has been selected, the negative reflection towards the m = -1 Floquet mode can be highly efficient (>0.9) over a relatively broad frequency range and for a broad range of incident angles. This was finally verified experimentally in a water tank for the rigid surface, using a metagrating build from cshaped brass tubes.

The broadband efficiency and simple design are significant steps towards real-world applica-tions in underwater acoustics, where most commonly encountered acoustic signals, from sonar to noise pollution, are broadband. Fabrication of acoustic metasurfaces is more challenging for water or water-like media than for air, since stiff and heavy materials are required to provide a large impedance contrast. We have demonstrated here that a simple brass hollow cylinder opened on the side is a viable structure for underwater applications. The broadband and broad-angle high efficiency in retro-reflection could be useful for telemetry applications based on pulse-echo methods.

The relation between the reflection angle and frequency could be used as a benefit to perform real time processing operations such as e.g. filtering, frequency discrimination (i.e. spectroscopy) or multiplexing / demultiplexing in underwater acoustic communications. The use as a noise barrier can also be envisaged in shallow water, where low frequency propagation is governed by guided modes reflecting back and forth between the water surface and the ocean bottom at grazing angles.

A metagrating, placed close to the water surface could have the ability to send back waves toward the source or toward the sea bottom with normal incidence, effectively blocking propagation. The design can be scaled to other frequencies, and applied to a free boundary. All analytical and numerical calculations were done under a linear acoustic propagation / linear elasticity assumption, which corresponds to the presented experimental conditions. Extremely high acoustic pressure levels may alter the behavior of the grating, by modifying the resonant behavior of the structure or through non-linear interactions of multiple scattered waves. This was not investigated here.

In naval applications, one usually seeks to reduce the acoustic power reflected back to the source by a target, which is precisely the opposite of what was achieved here. Further research is needed to develop other configurations, e.g. redirecting a normally incident wave toward a grazing direction. Using several atoms per period [START_REF] Torrent | Acoustic anomalous reflectors based on diffraction grating engineering[END_REF][START_REF] Fan | Metagratings for waterborne sound: Various functionalities enabled by an efficient inverse-design approach[END_REF] could potentially address this challenge. Finally, c-shaped atoms are known for providing strong Willis coupling [START_REF] Quan | Maximum Willis coupling in acoustic scatterers[END_REF][START_REF] Melnikov | Acoustic metaatom with experimentally verified maximum Willis coupling[END_REF], which was not discussed here but could also allow even more challenging wave control applications by canceling multiple Floquet modes simultaneously.

Appendix: Analytical model of wave diffraction by a periodic grating close to a reflective surface

We consider a periodic array of small particles whose center is at a distance h from a perfectly reflecting surface (either a rigid or a free surface). The configuration is depicted on Fig. 1(a). The surface is aligned with the x axis and at coordinate z = 0. The particles are at z = -h and the period is p. An acoustic plane wave of frequency f is incoming from the left, with and incidence θ i . The wave velocity and mass density of the propagation medium are c and ρ, respectively.

The incoming pressure wave is = ±p 0 e i(k x x-k z z) , (A.2)

p inc = p 0 e i(k x x+k z z) (A.1) with k x = k sin(θ i ), k z = k cos(θ i ),
where the (+) sign is for the rigid case and the (-) sign for the free boundary. Under a monopolar approximation (valid for small particles), the field scattered by the infinite array of particles is

p rigid/ f ree scat = A 0 +∞ ∑ n=-∞ H 1 0 k (x -np) 2 + (z + h) 2 ±H 1 0 k (x -np) 2 + (z -h) 2 e inφ , (A.3)
with H 1 0 the zero order Hankel function of the first kind and φ = 2π p sin(θ i )/λ a phase shift term related to the oblique incidence. The first term in the infinite sum is for the actual particles and the second term is for the "virtual" particles located on the other side of the perfectly reflecting surface. The Poisson's summation formula is then used to transform the above sum over an infinite number of particles to a sum over an infinite number of Floquet modes [START_REF] Tretyakov | Analytical modeling in applied electromagnetics[END_REF]. We obtain for the rigid boundary

p rigid scat = 4A 0 p +∞ ∑ m=-∞ exp i 2πm+φ p x-k zm z cos(k zm h) k zm , for z < -h, (A.4)
and for the free boundary

p f ree scat = -4iA 0 p +∞ ∑ m=-∞ exp i 2πm+φ p x-k zm z sin(k zm h) k zm , for z < -h, (A.5) with k zm = k 2 -2πm+φ p 2
the horizontal wave number for the m-th mode. The Floquet modes in the sums above are propagating if k zm is real and are evanescent if k zm is complex. Propagating modes leave the surface with an angle θ m r given by the relation k sin(θ m r ) = φ +2πm p , from which we can arrive at the grating equation sin(θ m r )sin(θ i ) = m λ p . To achieve perfect anomalous reflection, i.e. having the incident wave redirected with unit efficiency into a single plane wave whose angle can be controlled, we need two propagating Floquet and for the free boundary pinc + pre f = p 0 e ik z hp 0 e -ik z h = -2ip 0 sin(k z h).

(A.12)

Inserting these equations in eq. (A.10) and satisfying the conditions eqs. (A.8 -A.9) we find In both cases, we see that α must be purely imaginary for the m = 0 Floquet mode to cancel the specular reflection, which happens at the grating resonance [START_REF] Ra'di | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Tretyakov | Analytical modeling in applied electromagnetics[END_REF].

In addition, it is also required for perfect anomalous reflection that the power carried away from the surface by m = -1 mode matches the incoming power. The incoming power is equal to

P i = p 2 0 p 2ρc cosθ i . (A.15)
The power carried away by the m = -1 mode is, now simplifying notation using θ -1 r = θ r , for the rigid boundary One of these two last conditions, depending on the nature of the reflective surface behind the grating, has to be satisfied along the grating resonance condition in order to achieve perfect negative reflection. We can see that retro-reflection, θ r = -θ i is always a solution, for any distance h. However, as demonstrated in the main text, there are other non-trivial solutions and there exist values of h for which the conditions can be approximately satisfied over a broad range of frequencies and incidence angles, leading to a robust broadband negative reflection effect.

FIG. 2 .

 2 FIG. 2. Tuning the distance h between the grating and the rigid (a)-(b) or free (c)-(d) reflective surface.

  section IV below). Wave velocity and mass density of water are 1470 m/s and 1000 kg/m 3 , respectively. Modeling was done in the frequency domain. A single periodic cell of the grating was modeled as a rectangular box, using Floquet periodic boundary conditions to obtain the behavior of an infinite grating. A rigid (normal velocity = 0) or free (acoustic pressure = 0) boundary condition was set behind the atoms and a matched boundary condition at the opposite boundary was set to perfectly absorb the m = 0 and m = -1 reflected modes. The incident wave was a planar wave background pressure field with unit amplitude. The power radiated by each of the possible propagation mode was evaluated by integration over the opposite boundary, and expressed as(quadratic) reflection coefficient |R m | 2 .Results of the FEM modeling of the periodic grating for the rigid boundary case for h = 8 mm are presented on Fig.3. Panel (a) displays the reflection coefficient of the anomalous reflection m = -1 for all incidence angles from 0 to 90 • . Efficiency remains excellent (> 0.9) for a broad range of incident angles (14 • to 54 • ) around the design angle of 30 • . The bandwidth is shrinking for extreme angles but is still reasonably large angles from about 20 • to 40 • . The reflection angle here depends on both incidence angle and frequency, according to Eq. (1). Retro-reflection (θ r = -θ i ) is achieved along the purple line. More than 50% of the incoming power can be retro-reflected for frequencies ranging from 65 to 111 kHz, and angles from 23 • to 43 • . The sharp bright lines in panel (a) correspond to the appearance of new possible reflection modes among m = -2, m = -1, and m = 1. The reflection coefficients of each channel are plotted on panel (b) for θ i = 30 • . Until 115 kHz only the m = 0 and m = -1 modes exist, in agreement with Eq. 1. Above 70 kHz, the m = -1 reflection mode becomes dominant, up to 108 kHz. Between 74 kHz and 103 kHz, more than 90% power is directed toward the anomalous reflection direction, demonstrating near unit efficiency over a large bandwidth (≈ 1/3 of the central frequency). Panels (c) and (d) show respectively the incident and reflected wave field (real part of the acoustic pressure) at 87 kHz (p = λ ) for θ i = 30 • . In this case retro-reflection is perfect, with efficiency > 0.998. At 75 kHz and 100 kHz (Fig. 3(e)-(f)) anomalous reflection is less than perfect, as visible from the slightly distorted wave fronts, but remains efficient.

30 FIG. 3 .

 303 FIG. 3. FEM results for the proposed metagrating over a rigid surface with height h = 8 mm. (a) Reflection coefficient of the anomalous reflection (m = -1) channel as a function of incident angle and frequency. The purple line indicates retro-reflection (θ r = -θ i ). (b) Reflection coefficient of the channels m = 0 (specular), m = -1, and m = -2 for θ i = 30 • . The bandwidth over which |R -1 | > 0.9 is highlighted in grey. (c-f) Maps of the incident (c) or reflected (d-f) acoustic pressure field. At 87 kHz, retro-reflection is perfect. At 75 kHz and 100 kHz anomalous reflection is less than perfect but still efficient.

FIG. 4 .FIG. 5 .

 45 FIG. 4. FEM results for an acoustic free surface, with height h = 13 mm (a-b) and h = 2.5 mm (c-d). Maps (a) and (c) show the reflection coefficient of the anomalous reflection (m = -1) channel as a function of incident angle and frequency. The purple line indicates retro-reflection (θ r = -θ i ). Plots (b) and (d) show the reflection coefficient of the channels m = 0 (specular), m = -1, and m = -2 for θ i = 30 • . The bandwidth over which |R -1 | > 0.9 is highlighted in grey.

Fig. 7

 7 Fig. 7 provides an alternative representation of the experimentally recorded pressure at 87 kHz, at which the metagrating is the most effective, for the same four cases : (a) no metagrating, (b)-(d)

FIG. 6 .

 6 FIG. 6. Experimental results for the constructed metagrating. (a) Amplitude of the reflected pressure wave without grating (calibration). (b-d) Acoustic pressure reflected by the grating for three incidence angle (15 • , 30 • , 45 • ). (e-g) Profiles along the specular direction m = 0, and comparison to FEM results. For θ i = 30 • , the m = 0 mode is reduced by more than 10 dB over a broad frequency range (≈ 30kHz). The surface is still effective for θ i = 15 • and 45 • .

FIG. 7 .

 7 FIG. 7. Polar plots of the reflected pressure at 87 kHz, at which the grating is the most effective. (a) Without the metagrating and (b)-(d) with the metagrating for incidence angles 15 • , 30 • , and 45 • , respectively. The direction of the incident wave is displayed as a black arrow. Amplitudes are normalized to the amplitude of the specular reflection at 30 • .

  and k = ω/c = 2π f /c = 2π/λ . The pressure wave reflected at the boundary (in the absence of grating) is p rigid/ f ree re f

1 α = 2i p ωc cos θ i cos 2 (cos θ i sin 2 (

 122 kh cos θ i ), for the rigid boundary and (kh cos θ i ), for the free boundary.(A.14)

2 (

 2 kh cos θ r ) k 2 cos 2 (θ r ) .(A.17)Using the conditions eqs. (A.8-A.9) for p 0 above and equating incident and carried away powers lead to1 cos θ i cos 2 (kh cos θ i ) = 1 cos θ r cos 2 (kh cos θ r ),for the rigid boundary and (A.18)1 cos θ i sin 2 (kh cos θ i ) = 1 cos θ r sin 2 (kh cos θ r ),for the free boundary. (A.[START_REF] Epstein | Unveiling the properties of metagratings via a detailed analytical model for synthesis and analysis[END_REF]

ACKNOWLEDGMENTS

The authors thank H. Cahingt for his help with the maintenance and preparation of the experimental facility, and F. Duhamel for the preparation of the grating. This work received financial support from the Region Normandie (project number R2019-RIN-0027).

modes : one being the controlled reflected wave and the other one canceling the specular reflection. One of these two modes is necessarily the mode m = 0, as k z0 = k z is real. The second can be either m = 1 or m = -1 depending on the sign of θ i . We assume here that θ i > 0, in which case we need the modes m = 0 and m = -1 to be propagating and all other to be evanescent. A mode m is evanescent if sin(θ m r )sin(θ i ) = m λ p has no solution. For m = 1 no solution exists if p < λ 1-sin(θ i ) . Satisfying this condition also ensure that modes m = 2, 3, ... are evanescent. For m = -2 to be evanescent (and also m = -3, -4, ...), we get to the condition p < 2λ 1+sin(θ i ) . For p = λ , both conditions are satisfied for all angles θ i > 0.

The total wave field is then p tot = p inc + p re f + p m=0 + p m=-1 , which leads to, for the rigid boundary,

where k sin(θ

. For the free boundary we arrive at

For the m = 0 mode to cancel the specular reflection it is therefore required that

, for the rigid boundary and that (A.8)

, for the free boundary. (A.9)

The 0-order scattering coefficient A 0 is related to the effective monopolar polarizability α of the particle by the relation [START_REF] Melnikov | Acoustic metaatom with experimentally verified maximum Willis coupling[END_REF] A

where pinc and pre f denote respectively the incident and reflected pressure at the center of the particle. For the rigid boundary we have, for a particle located at x = 0 and z = -h pinc + pre f = p 0 e ik z h + p 0 e -ik z h = 2p 0 cos(k z h), (A.11)


Adrien Halnaut

Romain Giot

Romain Bourqui

David Auber

Compact visualization of DNN classification performances for interpretation and improvement

come

i "output" -2021/11/19 -8:22 -page 20 -#34 i i i

A. Introduction

Deep-learning [LBH15] based approaches are used in various contexts and dominate most historical methods, especially for classification problems. Even when training datasets are not large enough to train Deep Neural Networks (DNN), it is possible to use transfer learning with a pre-trained DNN by fine-tuning [GYA17] it, or by extracting features and feeding them to a conventional classifier [Tan13]. DNN can be represented by graph linking computational blocks between each other: each block processes the output of one or several previous one. Theses blocks are applying simple function defined during the model construction and its parameters are defined using weights. Such weights are data-dependent and computed during the training phase of the model. The black-box feeling is one of the largest issues. Indeed, even if each block is individually well understood mathematically, its behavior depends mainly on the training data (i.e., their impact on the learned weights). As a consequence, one can hardly know what processing these blocks are doing and why. However, it is well admitted that the model's first layers extract low-level features, while the latest ones extract high-level features specific to the application problem [EBCV09]. Two non-exclusive strategies can help to open this black box. (i) Explainable deep-learning where the architecture of the DNN emphasizes its explainability [ZNWZ18], even if it could negatively impact its performance, and (ii) Interpretable deep-learning where additional processes extract information by computing a more explainable model [ZYMW19], computing some saliency information [BML + 16, AFMB + 20] or using information visualization techniques [HKPC18]. Understanding, even partially, the classification process of a trained network helps in both understanding which data classes are easier to classify than other and where the network can be further improved for its task, in a different way than what layer pruning techniques do [HZS17].

This chapter is an extension of the paper [HGBA21], presenting a new method for interpretable deep-learning based on information visualization techniques. The task to solve corresponds to the analysis of the data classification over layers. It aims at analyzing how all input samples of a dataset are globally treated by any part of the network. In opposite to most papers of the literature on attribution-based methods, the focus is not for a specific sample but for a full testing dataset. The method allows focusing on successive layers that better (or worst) discriminate the samples. It also displays the complete network and the data behavior for each of its computational blocks.

It originality relies on the fact it focuses on both all samples and full architecture and has the advantages of (i) using less screen space than existing methods despite the amount of information to display; (ii) fitting to any network that can be represented as a directed acyclic graph; (iii) using the same encoding for input data, inner blocks and final result.

The remaining of the paper is organized as follows. Section B presents related works in visualization for DNN, space filling curves. Section C describes the proposed method. Section D provides the details of the experimental protocol. Section E discusses the results and provide directions for future work. We finally draw conclusions in Section F.

B. Previous Works

Our proposed method aims at visually interpreting how DNN behave using a spaceefficient method. For this reason, this section firstly presents previous works on deep neural networks visualization then focus on dense, pixel oriented visualization methods.

B.1. Visualization for the Interpretation of Deep Neural Networks

Some works focus on single views that can be reused in other works or embedded in more complex applications. GradCam [SCD + 17] aims at generating a heatmap for a single input to highlight the spatial location that greatly supports the final decision. It is computed thanks to the gradients from the logit of a target class up to the latest convolutional layer. When the input feature is an image, the heatmap can be straightforwardly visualized and understood on it. Other methods rely on different concepts to achieve the same objective such as LRP [BML + 16] on the concept of relevance, or other work [AFMB + 20] that only uses information collected during the forward pass. Instead of focusing on a single input sample, it is also possible to focus on the complete dataset. Some use Sankey-diagram analogy [HGBA20] to highlight the processing flows. Others project activations obtained at a specific layer in a 2d space to verify how the network sees the data at this specific point [RFFT17]. Such an approach is also common in the literature using T-SNE projection In opposite, several complete tools treat industrial problems. Some of them are generalist enough to be used in almost any scenario, such as Activis [KAKC18] (that focuses on the visualization and comparison of activation of a single selected layer), while some others are restricted to some specific networks or evaluation scenarios. CNNVIS [LSL + 17] is tailored for CNN and uses a visualization that relies on aggregation of layers (not all layers are depicted), filters (filters that behave similarly are grouped) and data (a subset of the samples are depicted). DQNVIZ [WGSY19] has been designed for Deep Q-Network explanation in the specific context of Atari play.

B.2. Hilbert Curve in Information Visualization

Dense pixel-oriented methods aim at improving both the data-ink ratio ?? and the visualization size by displaying a unit of information on a single pixel while avoiding unused pixels. Keim reviewed various pixel-oriented visualizations [Kei00] and asserts that space-filling curves, such as the Hilbert one [Hil35], are among the bests to project ordered elements in a screen space while preserving the distance of the onedimensional ordering in the two-dimensional arrangement. Blanchard et al. [BHL05] have shown that to display images, reduced to a one pixel representation, on an Hilbert curve produces coherent and identifiable clusters. Auber et al.

[ANM07] have also shown the interest of such visualization, when complemented by tailored interaction techniques, to explore datacubes of several dozen of millions of elements. Since these previous successes, we have selected the Hilbert curve to project our data in a square; a curve of order n contains 4 n elements [ANM07].

C. Proposed Method

Figure 3.1 describes the proposal with the "Nested blocks and guidelines model"[MSQM15] among various description levels: domain (who is concerned by which problem), abstraction (which data is used or generated to solve which task), technique (which methods are used) and algorithms (how these methods are implemented).

Additionally, Figure 3.2 lists the successive steps involved in the method. The requirements of the proposed method are: to be space efficient (R1) while displaying information from all samples (R2) in all layers (R3) of the network to solve the task "classification quality analysis over layers" (R4).

C.1. Domain Level

The proposed method fits the needs of networks designers and trainers that want to verify how the data is grouped by the various layers of their classification network. From the analysis of these groupings, they could infer hypotheses that aim at being verified with other techniques. Such hypotheses are related to input sample properties and network errors. Non experts would better understand how DNN work by looking at the representation of simple networks and datasets.

C.2. Abstraction Level

The proposed method considers an already trained DNN N with a compatible test dataset D test . N is a network (i.e., graph) of operations (i.e., nodes) N = (O, E). Its sources s • ∈ O are the identity function on data input (i.e., samples) and its sinks t • ∈ O are its outputs (i.e., classes probability). N has multiple sources for a multimodal system, but always a single sink as we are restricted the use case of standard classification. The other nodes o • ∈ O \ {s • ∪ t • } correspond to any operations (e.g., Each sample d i ∈ D test is fed into the network and the output (i.e., activations) of each operation o j is stored in a j i ; we assume operations are ordered depending on the execution flow. These activations consist of tensors whose order depends on the underlying operation and whose dimensions size depends on the input data of the network. Each operation o j consumes at least one result a k i |k < j computed by a previous operation except for the sources where a • i corresponds to the raw data (of the targeted modality in a multimodal scenario). Thus, a sample d i is represented by a set of activations A i = ∪ j {a j i } and the complete dataset D test is represented by an ensemble of sets of activations D activations = ∪ i {A i }.

The activations can be optionally preprocessed to fall within compatible domains as their domain is not controlled: D preprocess = i j Preprocessed a j i . This preprocessing method is a parameter of the workflow.

C.3. Technique Level

As the method aims to display (a) the dataset and its groundtruth (R2), (b) the architecture of the network (R3) and (c) its impact on the complete dataset (R4), we propose an encoding relying on both D preprocess and N.

Groundtruth encoding. The groundtruth of the dataset is depicted with a legend where each class is represented by a colored rectangle (sample encoding) followed by a black text (classe name).

Network Encoding. It is straightforward to layout N operations with a graphdrawing algorithm tailored for Directed Acyclic Graphs (networks are always DAG). Such technique is common in the literature [KAKC18, WSW + 17] and aims at computing the coordinates of each node (operation) in a plane while emphasizing the order of operations in the computing flow. Each node is depicted by a glyph that represents the whole dataset as viewed by the network at this specific operation. Thus, a specific encoding is used to map the activations i {Preprocessed(a j i )} of each node o j in the screen space.

Similarly to Ganlab [KTC + 18b], a dotted line is drawn between nodes that represents consecutive operations; the flow of data is revealed by the dots moving in flow direction. Some networks can be very deep with successive layers that do not bring additional information because they consist of data reordering. We allow the user to request the visualization of a simplified network where the corresponding nodes are removed (thus, their successors are linked to their predecessors), as such information brings noise to the representation. No special encoding is used to represent this information shrinking.

Samples encoding. As mentioned, we have chosen a pixel-oriented technique that relies on fractal curves (R1). For a given node o j , a maximal linear arrangement method is used to order the representation a j i of each sample d i in such a way that samples are positioned closely in the ordering according to a distance function. We assume close samples in the output space of o j corresponds to samples treated similarly by the network (i.e., considered to be similar). Once the samples are ordered, they are projected into a discrete pixel grid using a fractal curve that respects proximity relations. This way, screen space usage can be maximized (1 pixel per sample) and we are assured that close samples are drawn closely on the screen (however close pixels on the screen are not necessarily close in data space). Two visual encodings can represent this curve. The first one, absolute coloring, explicitly draws samples of each class with the same color. The second one, relative coloring, uses a gray-scale to emphasize label difference between adjacent nodes and identify zones where different labels are present. It can be used de facto when the number of classes is too high to be efficiently discernable by a human using regular class colorization. When using the absolute coloring scheme, the user can choose to only visualize a specific class to observe the spread of its samples over the layer. The name of the layer is written above its fractal representation, and a quality metric (presented later in this chapter) is written below it.

C.4. Algorithms Level

The model topology is drawn using the well-known Sugiyama [STT81] algorithm and each node is depicted with a specific fractal-based glyph that represents the ordered samples. The Euclidean distance is used to compare the activations generated for all the samples on the same operation. It reflects the dissimilarity between samples in the Euclidean space; we consider that each neuron activation has the same impact as others in the full network processing. These distances are then compiled into a n × n sized distance-matrix, n being the number of compared samples. In real use case, neurons have different impact on the final prediction than others. Some pre-or post-processing methods, such as the LRP [BML + 16] method as done in [HGBA20], can be applied to the activation maps in order to reflect that behavior. However, we decided not to apply those methods because of the unsure interpretation on model topologies using branches, such as our chimeric DoubleLeNet5 (section ??) or the widely used ResNet [HZRS16] which use residual connections. Using ordering methods [BBHR + 16], data can be ordered in a queue with similar elements placed next to each other using their dissimilarity matrix. By using the VAT algorithm [dSW18] on the dissimilarity matrices, we found a progressive definition of similarly processed samples, resulting in clusters (or "black squares" as defined in the original paper [dSW18]) reflecting the progressive recognition by the model over the layers we attempt to show. The order computed by this algorithm can then be applied on a 1d-space to display similar data indexes next to each other. Using a fractal curve, we transformed this 1d-space into a 2d-space which is more suitable for data visualization. The fractal curve chosen to map each sample into a pixel-grid is the Hilbert curve [Hil35] because of its ability to place points in a discrete space (this is not the case of Gosper curve [Gar76]) and the absence of "jumps" in the curve (this is not the case of the Z-order curve [Mor66]) which ensures that two consecutive samples are adjacent. The order in which each sample is positioned is following the same order computed by VAT on the previous step. When the number of test samples is lower than the number of pixels available in the curve, we skip half of the missing positions in the beginning of the curve (and thus half of the missing positions at the end of the curve); making a "hole" in the curve but keep the sample centered in the glyph.

In the absolute coloring, each pixel sample is being colored according to its groundtruth class, which is different for each class. In the relative coloring, the colors depend on the number of similar labels for the pixel of interest in its sample ordering. That gives three possible values (0 for an outlier with no neighbors of the same class, 1 for a previous or next label different, and 2 when the three successive samples are of the same class). The absolute colors come from a palette of diverging colors while the relative colors or black (0), gray (1) and white (2). Computing in the ordering space instead of the picture allows to no highlight the visual border inherent to the fractal curve. Placing the cursor on a class in the legend selects this specific class and draws only its samples with the appropriate absolute color.

The machine learning community provides various evaluation metrics (e.g., accuracy or cross-entropy) to evaluate the quality of the network by comparing its output to a ground-truth. By definition, they cannot be applied at each layer, but we still need to provide hints to the user of their efficiency. We have defined a quality metric, based on the local homogeneity of the layer's visual representations, which counts the number of neighbors of a given pixel that are of the same color (i.e., the number of samples that belong to the same class). We normalized it between 0 and 1 to ease its comparison (however, as the normalization does not consider the mandatory borders, 1 is an unreachable value). We assume that to quantify the quality of the visualization is strongly related to the ability of the layer to separate data.

D. Experimental Protocol

Several scenarios, that rely on a test dataset and a trained network, illustrate the efficiency of the proposed method. In this section, we present datasets and networks that we used in our evaluation.

Datasets. MNIST [LeC98] is a standard dataset used in handwritten recognition from 28 × 28 grayscale images. Even simple networks are able to perform almost perfectly on this 10-classes dataset. We use it to illustrate classification on easy data. Fashion-MNIST [XRV17] shares a similar distribution as MNIST but is composed of images of clothes instead of digits. Classification performance is usually lower than with MNIST. We use it to illustrate classification on averagely difficulty dataset, closer to actual classification problems.

Both datasets are composed of 60 000 samples to train the model and 10 000 samples to evaluate the model.

Networks. In our evaluation, we make use of three networks summarized in Table 3 • LeNet5 [LBBH98] is a simple and historical CNN that provides good accuracy results on MNIST. Its topology is simple enough to get a grasp on how data is being transformed across the model. It is also easy to train with its low parameter count, but that simplicity comes at the cost of lower accuracy results in more complex recognition tasks. • DoubleLeNet5 is a chimeric network we have created to illustrate the ability of the system to handle networks with several branches. It consists of two LeNet5 minus the prediction layer that process in parallel the same input data, but one of the branches input has an image rotation step applied before being processed by the convolutional layers. The two branches are then concatenated before being fed to the prediction layer. The image rotation step is not represented in the DAGs has it has not be implemented in the same way as other layers during the evaluation.

Performance wise, this model targets the same kinds of data as LeNet5, with a minor performance gain. • VGG16 [SZ14] is a deep CNN usually used on complex datasets composed of large color images, with a thousand recognizable classes, such as ImageNet [RDS + 15]. Its robustness allows it to reach fairly good accuracy results on target tasks, but comes with a heavy computation cost and cannot be trained in a reasonable amount of time on standard computers. In this paper, the convolutional blocks of the VGG16 model are already pre-trained with the ImageNet dataset, and are not retrained when training the prediction layers.

D.1. Couples of Network and Dataset

We have selected meaningful combinations of network and dataset.

• Easy scenario: LeNet5 uses MNIST which illustrates a well performing system.

• Generalization scenario: LeNet5 uses Fashion-MNIST which illustrates a system with more classification errors. • Branch scenario: DoubleLeNet5 predicts Fashion-MNIST which illustrates a usage case with non-linear network architecture. • Simplification scenario: VGG16 is processing MNIST. It illustrates the use of a complex network to solve a simple task. By applying the visualization pipeline and observing resulting glyph, we propose a straightforward improvement for the model in both accuracy and complexity for the model.

D.2. Implementation and Execution Infrastructure

The TensorFlow framework [AAB + 15] is solicited along with Keras. [C + 15] to train the studied models with said datasets. Each layer output, processed as potentially very large high-dimensional data, are saved into machine cluster handling the dissimilarity matrix computation, which makes use of a large pool of memory (around 2 Terabytes in our infrastructure). The resulting matrices are small enough (for our experiments) to fit and be processed on a recent laptop. The matrix manipulations of the original VAT algorithm [dSW18] are implemented using the ArrayFire library [YAM + 15] for their efficient matrix computation abilities. This part of the process only produces the data for the visualization tool and thus can be seen as a backend infrastructure. Fractal images (i.e., glyphs) are then generated by relying on the Rusthilbert library [Ski04]. The visual and interactive part corresponds to an HTML application written in Typescript relying on D3.js for the visualization, D3-dag for the Sugiyama implementation and webpack for the build system.

E. Results and discussion

The complete results are accessible at the following address: https://pivert. labri.fr/frac/index.html. We strongly recommend to view the results online as the images here are severely undersized on the paper where 1 pixel represents several samples) Fig 3 .3 depicts still resized representations of the proposed method for several scenarios, while their confusion matrices are presented in 3.4.

E.1. Results

The accuracy for the system of the easy scenario is 98.96%. This is clearly reflected by the heterogeneous organization of the 09-prediction glyph in the visualization. Looking at the successive operations that correspond to activation functions (01-relu0, 03-relu1, 06-relu2, 09-prediction), we observe an improvement in the quality of the representation, and thus better discriminability ability over layers.

Classes 1 and 0 seem to be discriminated early in the the network and can be considered as "easy" classes for the model. Several elements of the 6, 7 and 8 classes are also recognized early. The first dense layer 05-dense0 brings a dramatic improvement of the discriminability in the classification. The accuracy in the generalization scenario is 73.44%. Compared to the prediction layer's glyph from the easy scenario, this one is saltier. The previous layers overall are also less organized, which reflect the overall lower performances of the model on this dataset. The T-shirt and Trouser classes are differentiated early in the classification process while Ankle boots or Sandal are discriminated only at the end of the network.

The simplification scenario has the accuracy of 97% (lower than LeNet5 despite the more complex architecture). During the progressive classification, we can notice a succession of improvement and decrease in sample organization in the glyph. We assume this less efficient classification come from several reasons: they are pre-trained and not specialized for the task, and/or the model is too deep and layers are redundant. In the next subsection, we propose an improvement of the network architecture based on this observation. The realistic scenario illustrates the ability to draw networks with branches. The 73.75% accuracy of the model is very close to its linear counterpart in the generalization scenario. We observe the same tendencies in both branches as well as a similar sample organization than in the generalization scenario. They are also confident on the same classes.

E.2. Simplification Improvement

The glyph representing the fifth convolutional block of VGG16 17-block5-pool shows a large degradation of the classification process, with fewer subsets of data being similarly processed than at the end of fourth convolutional block 13-block4-pool. We imply that data is being over-processed by convolution operations which degrades the model performances. We propose two ways to improve the model architecture based on this assumption :

• VGG16-B4: Remove the fifth convolutional block altogether and directly connect the fourth block to the top of the model. (i.e., dense layers) • VGG16-B4+: Remove the fifth convolutional block but keep the Pooling layer before connecting the model's top. This ensure that the number of weight parameters for each neuron of the next layer is not increased after the modification. In both cases, the model's top is being re-trained from scratch (with only results from ImageNet processing for the convolutional blocks) to ensure that any noticeable improvement is due to the architecture modification instead of possible transfer learning from previous iteration.

The results of VGG16's simplification is shown in Table 3.2. Removing the fifth convolutional block of VGG16 improved the classification performances of the model while reducing its complexity. Which improvement is better for the use case is up to the user's decision. However, even without knowing the parameters of each VGG16's layers nor the nature of the processed data, the user can notice the classification degra-

E.3. Discussion

Compared to widely used t-SNE projection [RFFT17], the method presented in this chapter has a fairly more efficient use of the screen space. Furthermore, all of the layer activations can be displayed on the same screen space without overploting data. Focus on pixel-scale usage is emphasized, but not all points of the curves are used; black pixels correspond to unused pixels because test datasets are smaller than what is technically possible with such display size. Indeed, in the case of Hilbert's curve, only datasets of size 4 n can entirely fit into the curve. In our experiment, 4 7 -10000 pixels are lost, which is roughly 39% of the picture for each layer. It is thus possible to evaluate larger datasets without using more space on the screen, resulting in better data-ink usage. Another observation is the effectiveness of samples projection over the fractal curve to depict the classification performance over layers. Usually, representation of sample ordering is getting better over layers which means the network is progressively better at separating classes of samples. Dataset and subsets of dataset classification difficulty are also represented: Fashion-MNIST, which is a problem more difficult than MNIST, is thus less well organized. By construction, the very first node corresponds to the projection of the raw dataset; the noisier it is, the more complex it is to distinguish its samples without extracting additional features. The representation clearly depicts this point and its quality metric is worst for Fashion-MNIST than MNIST. The very last node corresponds to the projection of the softmax values; the noisier it is, the worst the network's accuracy is. The final representation is complementary of a confusion matrix (see Figure 3.4) as it provides more information about classification efficiency of the models.

A labeled dataset is currently needed to color the pixels. It limits the use of the method to a test dataset and not a real world unlabeld dataset. However, it is still Focusing on a specific class helps to track the evolution of the sample processing for that chosen class over layers. Figure 3.6 illustrates a the oversizing of the VGG16 network on MNIST by observing that samples of selected class tend to be processed similarly at the 09-block3-pool layer, whereas it is not the case anymore around layer 17-block5-pool.

E.4. Future Work

The method is resource consuming, mainly due to the need of storing dissimilarity matrices in memory, which is of size N 2 , with N being the size of the tested dataset. As a future work, it would be interesting to study whether one could use smaller part, or estimation, of the dissimilarities to approach a similar visualization. Additionally, the ordering of the samples highly depend on the Euclidean distance that is known to not be efficient in high dimensional spaces; other metrics need to be compared. The approach is satisfactory using interaction, but is not yet self-sufficient. Indeed, it provides a good overview of how the classification is handled but lacks of interactions to track the progression of a single sample or group of samples (in opposite to our previous work that specifically focus on this point [HGBA20]) in the network. Such investigations have to be held; for example, some sort of consistency in the sample position between two successive glyph would help in tracking elements. Furthermore, being able to focus a single sample instead of the whole class would help in determining the cause of miss-classification by the network (e.g. the model made confusion between a 6 and a 8 at the first convolutional layer, which led in miss-classification for the rest of the processing). The Hilbert curve is very efficient to place the samples in its reserved space. However, there is a high probability that the number of elements in the dataset to visualize is lower than what is possible with the curve. It would be interesting to implement additional interactions that use this additional space; or use grid-based projection methods instead of fractal ones. To subsample or sample with replacement the dataset with a number of samples equals to the curve length, and that follows data distribution, could also be interesting. The standard Sugiyama algorithm does not consider the screen space size; a modified method should be used in order to project the graph on the screen in a way that does not necessitate to horizontally scroll the screen to see it [LLS + 18].

F. Conclusion

Deep learning classifiers are progressively replacing handcrafted and understood standard classifiers for various fields. This significant gain in performance and accuracy is counterbalanced by a steep difficulty in understanding how and why they perform so well. Information visualization is one solution to this fill this lack of interpretability. We have presented a pipeline consuming a trained network and a dataset which produces an interactive representation depicting both the network's architecture and the behaviors of each layer when they process the test samples. Such system allows to visually analyze the classification quality over layers of a dataset and could be used to visually detect patterns in the data. This analysis would lead to a hypothesis about the performance of the network. However, such hypothesis would need then to be verified by other means.

This approach has been validated on various scenarios and shows its interest and limits that could be overcome in the future. Extensions with various specific interaction methods to also focus on individual data, efficient data subsampling and dense pixel-based glyph construction with better screen-space usage and/or less restrictions would imrpove the method for more complex and precise network analysis.

  Figure 3.1: Nested blocks and guidelines [MSQM15] representation of the proposed method. Dotted italic blocks corresponds to existing ones; plain straight blocks are defined in the work.

Figure 3

 3 Figure 3.2: Summary of the proposed method. Dataset D is fed to the network N. All activations are collected, eventually preprocessed, and finally ordered at each layer.The ordered samples are drawn along a fractal curve at each operation of the network that is placed on the screen using a graph drawing method.

  Compact Visualization of DNN Classification Performances for Interpretation and Improvement 25

  Compact Visualization of DNN Classification Performances for Interpretation and Improvement 27

  Compact Visualization of DNN Classification Performances for Interpretation and Improvement

  Figure 3.3: Illustration of results on some scenarios. The simplified version of the network is drawn. The confusion matrices are presented in 3.4 for comparison with a standard visual evaluation method. Larger images are available on the website https: //pivert.labri.fr/frac/.

  Figure 3.4: Confusion matrices of the systems presented in Figure 3.3.

  Figure 3.5: Comparison of the absolute and relative color schemes. No data is depicted in black for absolute and light blue for relative color schemes.

Figure 3

 3 Figure 3.6: Analysis of the spread of samples of class 2 over layers. Such representation could indicated an oversizing of the network by looking at the separation effect around layer 17-block5-pool.

  Methods and applicationsan important drawback remains: the representation is not space efficient and there is no guarantee that overlap does not occur. Our proposed method solves these two issues. Other works create applications for educational purpose to visually explain how some specific deep systems perform. For example, Tensorflow playground [SCS+ 17] focuses on simple DNN, CNN 101 [WTS+ 20] focuses on CNN, Ganlab [KTC + 18a] focuses on Gan system and Adversarial Playground [NQ17] illustrates the concept of adversarial examples. Even if they are visually appealing, these systems can hardly be used for industrial scenarios.

	i	i
		"output" -2021/11/19 -8:22 -page 22 -#36
	22	Explainable Deep Learning-
		[MH08]; however,
	i	i

Table 3 .

 3 1: Number of layers, parameters and activations per sample for each network.

	Network	Layers Parameters Memory Usage
	LeNet5	10	1 182 006	447.16 ± 19.7 MB
	DoubleLeNet5 18	1 646 370	54 570
	VGG16	18	33 638 218	308 244

ACKNOWLEDGMENTS

his work has been carried out with financial support from the French State. Experiments presented in this paper were carried out using the Labo's in the Sky with Data(LSD), the LaBRI data platform partially funded by Region Nouvelle Aquitaine.


Laurence Miègeville
email: laurence.miegeville@univ-nantes.fr

Patrick Guérin
email: patrick.guerin@univ-nantes.fr

Dominique Dreulle
email: dominique.dreulle@intradef.gouv.fr

Investigation of Correlations Between Driving Patterns and Power Demand of Auxiliary Devices Aboard Military Vehicles

Keywords: Auxiliary equipment, driving cycle, hybrid electric vehicle, Markov process, military vehicle. I

With the broad introduction of ever more electric equipment and devices aboard military vehicles, demands in energy supply is being increased to fit the mission requirements. The original aspect of this work is to propose a comprehensive review of the correlations between the energy consumption of the auxiliary equipment and the driving patterns in military areas. The tendency to serial correlation of the auxiliary power is also examined and argued. The longer term purpose is to determine the best suited model of the auxiliary power in order to properly size the electrical sources such as storage devices and generators.

INTRODUCTION

Following the example of the civil automotive sector, the electrification of military vehicles has been widely increased for the last fifteen years. A full range of auxiliary equipment is now electrically actuated or powered such as tank turrets, air conditioning systems, gear boxes, or cooling strip; also new equipment are currently introduced like optics detector (OD), laser warning receiver (LWR) or else identification friend or foe (IFF) devices. Then, the addition of hybrid-electric propulsion can be expected to require on-board electrical power distribution systems able to reach up to hundreds of kW. Besides in a distant future, electrical armaments should be installed onboard, which means excessively high aboard power level. In the medium term however, electric propulsion will remain the most energy-intensive.

The interest in promoting electric drives is not solely driven by lower fuel consumption or beneficial consequences for the environment which remain a key priority for the civil sector, but more especially for the operational advantages that could be derived from. A good example of this is the likely introduction of a stealth mode when the vehicle is only powered by the electric propulsion. Some other benefits are also expected, like the availability of over-torque or transitional overpower in order to climb over natural barriers, the reduction of the thermal signature when running on a all-electric propulsion mode, or else a greater flexibility in vehicle design following a removal of mechanical linkages. Finally, the full hybrid system should enable greater relevance to the increasing electric power consumption required by the auxiliary systems and new electric equipment aboard.

The purpose of this paper is to propose a comprehensive review of the correlations that are likely to exist between the power demand of the auxiliary devices onboard and the driving parameters for a military vehicle typically used to move forward in numerous and specific theatres of operation.

II. REVIEW OF THE DRIVING CYCLES USAGE

A. Military Usage

Some of the research dealing with the High-Mobility Multipurpose Wheeled Vehicle (HMMWV) of the US army offers relevant insights regarding the optimal power sizing of the energy sources. The approach and scope of those studies mainly focused on the implementation of behavioural multiphysical models for the hybrid powertrain.

A parametric approach implemented on the Advanced Vehicle Simulator (ADVISOR) was proposed in [START_REF] Antonis | Modeling and Simulation of Various Hybrid Electric Configurations of the High-Mobility Multipurpose Wheeled Vehicle (HMMWV)[END_REF] so as to compare some of the vehicle's performances between different hybrid motor powertrain concepts. The analysis is then conducted by varying the rate of hybridization defined by the ratio between the electric and thermal power demands. In [START_REF] Bargar | Modeling and Verification of Hybrid Electric HMMWV Performance[END_REF], a modelling approach applied to a series HEV, namely, the U.S. Military's HMMWV XM1124 was developed from the PNGV System Analysis Toolkit (PSAT) of Matlab/Simulink. Validated on the basis of the FUDS drive cycle, the model should be however complemented by tests with prototypes. Those studies highlight the needs for having access to more representative driving cycles that can match with the military operating conditions. For example, the modelling of the auxiliary equipment is carried out from default values available from PSAT.

Following the previous research studies, [START_REF] Ducusin | Modeling of a Series Hybrid Electric High Mobility Multipurpose Wheeled Vehicle" Vehicle Power and Propulsion[END_REF] shows the possibilities for implementing with a far better accuracy a multi-physical model for a series-hybrid drive of a HMMWV of the US army by using the Ansoft Simplorer software. For that purpose, the author refers to the Urban Dynanometer Driving Schedule (UDDS) combined with the SAE J227a Schedule C drive cycle as input data, also considering the crossing of grades up to 60%. The auxiliary equipment modelling is however not taken into account. Following the same topic, [START_REF] Fish | Modeling Hybrid Electric HMMWV Power System Performance[END_REF] proposes the model of two hybrid drives for HMMWV prototypes from a specific tool developed by the Defence Advanced Research Projects Agency (DARPA). The studies aim at validating the basic models of the components available from that toolbox used for modelling the powertrain architectures of future vehicles. Then, basic tests on acceleration and grade-ability are performed and a comparison between simulations and experimental results is analysed on the basis of experience acquired through the testing of prototypes. Finally, a novel modelling and validation approach for hybrid electric vehicles (HEVs) proposed by [START_REF] Bogosyan | A Novel Model Validation and Estimation Approach for Hybrid Serial Electric Vehicles[END_REF] offers an improved simulation accuracy over most map-based modelling techniques. Applied to the US military's HMMWV XM1124, the results motivate the use of the approach for the estimation of variables unavailable for measurement in a variety of systems, such as engine torque, battery state of charge, generator torque, motor torque or fuel consumption.

B. Civil usage

Within civil applications, standard driving cycles are widely used for simulating the operating conditions of hybrid drives based on various sizing approaches to both reduce fuel consumption and decrease gas emissions.

In this respect, [START_REF] Desai | Optimal Design of a Parallel Hybrid Electric Vehicle using Multi-Objective Genetic Algorithms[END_REF] proposes a multi-objective optimization problem formulated for a parallel HEV, considering powertrain components, fuel economy and emissions minimisation, as design objectives. Performance of the designed parallel HEV is then evaluated using ADVISOR and considering the UDDS as a test case, in order to propose solutions for future designs. In [START_REF] Buecherl | Verification of the Optimum Hybridization Factor as Design Parameter of Hybrid Electric Vehicles[END_REF], the hybridization factor is advanced to an optimisation parameter of the required electric power at the stage of design for a hybrid-parallel propulsion. The vehicle behaviour is then simulated considering both the Federal Test Procedure FTP-75 and the New European driving cycle (NEDC). The NEDC cycle is also applied by [START_REF] Wang | Experimental Characterization of a Supercapacitor-Based Electrical Torque-Boost System for Downsized ICE Vehicles[END_REF] to perform a downsizing strategy by reducing the energy consumption of an internal combustion engine combined with an electric motor powered by a super capacitor.

Reference [START_REF] Wang | Design and Simulation for a Series-parallel Hybrid Electric City-bus[END_REF] uses a logic-threshold power management strategy based on the optimal engine torque for a Seriesparallel Hybrid Electric city-bus powertrain. The simulations are then performed on the basis of the ECE drive cycle, and component selecting and parameter matching of the powertrain system are implemented for reduced fuel consumption purposes.

All of those approaches refer to standard driving cycles for being tested and validated. But using such schedules for the military area shows two major drawbacks: firstly, the existing drive cycles are mainly representative of a civil application and secondly they typically cover a single vehicle usage. In order to provide proper answers to this lack of representativeness also noticed for civil applications, a number of studies have chosen to promote stochastic approaches able to generate a large number of driving patterns for covering a larger variability of the vehicle usage.

In this respect, [START_REF] Ripaccioli | A Stochastic Model Predictive Control Approach for Series Hybrid Electric Vehicle Power Management[END_REF] proposes a stochastic approach dealing with power management for a series-hybrid drive. The power required by the vehicle in motion is then assessed by a Markov chain which matrix of transition probabilities is implemented from standard driving cycles representative for urban, suburban and expressway driving conditions. A Markov process is also considered in [START_REF] Lin | A Stochastic Control Strategy for Hybrid Electric Vehicles[END_REF] so as to generate power demand patterns and to assess the speed of rotation of the wheels for trucks using a parallel hybrid power system. The generated statistical cycles are then used so as to implement a comprehensive energy management strategy between the both power sources combining a heat engine with a battery to provide energy. The optimisation process is then applied to maximise fuel economy and minimise the environmental impact with a subsequent reduction of the greenhouse gas emissions. Some other works still dealing with energy management prefer a dynamic programming method [START_REF] Kolmanovskyt | Optimization of Powertrain Operating Policy for Feasibility Assessment and Calibration: Stochastic Dynamic Programming Approach[END_REF] using a model based on a Markov process and applied to the NEDC driving cycle. Nevertheless, by still working on standard drive cycles, the major drawback is to limit the number of states in the Markov matrix and to reduce consequently the representativeness of the randomly generated cycles. To get round this problem, [START_REF] Johannesson | Assessing the Potential of Predictive Control for Hybrid Vehicle Powertrains Using Stochastic Dynamic Programming[END_REF] proposes a Markov process defined by a large amount of measured data including the vehicle speed and its location. More recently, [START_REF] Souffran | Simulation of Real-Word Vehicle Missions Using a Stochastic Markov Model for Optimal Powertrain Sizing[END_REF] promotes a new methodology for modelling real-world driving patterns by assigning to the state as defined in the Markov formulation, a combination of three dependent parameters representative of the vehicle driving conditions i.e. speed, acceleration and road gradient while preserving the correlation between the three variables.

Apart the work of [START_REF] Han | Characterization of driving pattern and development of a driving cycle in a military area[END_REF] that proposes a multivariate statistical analysis to develop a gradient-sensitive driving cycle for military vehicles in operation, the use of the previous methodologies based on stochastic approaches is not yet widely applied within a military context. Also, the existing works are more concerned with applications to powertrain sizing and have consequently no real interest in the contribution of the auxiliary consumption to the energy requirements over a mission. Thus, the idea of this paper is to consider the role of the auxiliary equipment in the power demand of a military vehicle in accordance with its driving cycle. The authors propose then to investigate the correlations that are likely to exist between the variables associated with the driving patterns and the energy consumption of the aboard auxiliary equipment.

III. MAJOR FEATURES OF A MILITARY VEHICLE MISSION

Over its life cycle, a combat armoured vehicle is mainly devoted to tactical combats in conflict areas during external operations or for training exercises. Fig. 1 displays the driving pattern of an armoured vehicle over one full mission of a training exercise. It is then observed that the graph combines a number of stop and running sequences of quite different durations. The stop phases typically match with quiet mode or monitoring activities that are strong consumers of energy for supplying the aboard auxiliary equipment. These remarks are supported by the values of Table I. which illustrates the main characteristics of four successive missions of a military vehicle.

It can be then noticed some very specific features compared to a civil use. Firstly, the distance travelled in military areas remain quite small; secondly, the maximal velocity achieved over the mission appears far lower than the levels reached by a civil car. Also, a quite low average speed can be observed due to frequent stops with long durations. Varying between 2 s to 10500 s, the stop phases are indeed about 80% of the time of the mission. In this regard, the power demand during the stop phases of the vehicle in start and stop or stealth mode can no longer be overlooked, given the above-mentioned features of the vehicle stops. Finally, the influence of the road gradient already stressed in [START_REF] Souffran | Simulation of Real-Word Vehicle Missions Using a Stochastic Markov Model for Optimal Powertrain Sizing[END_REF][17] is also noticeable with sharp levels ranging from -20.76% to 17.94%. IV. CORRELATION ANALYSIS

A. Introduction

The set of equipment and auxiliaries available in a military vehicle include all electrical energy-using devices which do not belong to the propulsion chain. A number of examples can be mentioned, i.e. deployable turrets actuators, observation and sighting devices (cameras, light intensification optics, ..), user interfaces, onboard computers, door actuators, communications systems … This class of equipment contains more than 120 main consumers. Given the complexity of their functioning and the subsequent difficulty in considering a suitable behavioural model, the power demand will be modelled by a stochastic approach.

It is reasonable to expect that the energy consumption of some devices is correlated with the vehicle motion and the geometric pattern of the route. In practice, when the fencer is following a target with a sighting system, the turrets movement is clearly dependent on the trajectory, the distance and the target's travel speed. However, it may also depend on the path and speed of the vehicle in motion. The turrets actuators also are likely to be stressed according to the vehicle's own moving. Another example concerns the stabilizing mechanism for a weapon in shooting operation when the vehicle is moving: correlated with the vehicle speed, the surface roughness of the terrain is expected to generate dynamic effects on the weapon that should be controlled. When considering the time waveforms of the auxiliary power and the vehicle speed in Fig. 2 and Fig. 3, this assumption seems to be valid. It must be however confirmed on the basis of a more objective and stringent mathematical criteria.

The analysis of the correlation between the time-dependent variables is conducted using the Pearson Product Moment Correlation together with the Spearman's rank correlation. The scatterplot typically used for a first stage analysis provides a visual distribution of the data and tells about the nature of the linkage between two sets of variable.

B. Inter-correlation analysis

The physical quantities to examine are the variables linked to the vehicle mission like the speed, the acceleration and the road gradient as well as the power demand of the auxiliary equipment, respectively referred to as v, γ, α and P aux . The idea is to perform a sensitivity analysis to decide whether it is appropriate or not to correlate the use of the auxiliary consumers and the main mission on the ground.

Signs of correlation between the vehicle speed and the electric consumption of the auxiliary equipment can be directly detected from the recordings of Fig. 2 and3. The apparent relationship observed over the full mission seems to be obviously associated with the vehicle stop phases, which is supported by the calculation of the relative correlation coefficient, as shown by the comparison in table II between two configurations: the full mission (including stop phases) and the mission without stop sequences.

1) Vehicle in motion

By omitting the periods when the vehicle is stopped, correlation between speed and auxiliary power consumption is setoff, as confirmed in the scatterplots of Fig. 4. Conversely, as outlined by the correlation coefficients of Table III, no consistent evidence of linear and non linear correlation is observed between the auxiliary power demand and the other mission parameters that are the vehicle acceleration and the road gradient.

2) Vehicle stopped Besides, over stop sequences, the auxiliary power pattern looks like different according to the state 'on' / 'off' of the engine, as shown in Fig. 5. At zero speed with the engine switched off, the vehicle goes into a generally long-term sleep mode, where the use of the auxiliary equipment is quite different (small displacement of the weapon system, a little moving observation …). At zero speed with the engine running, the vehicle is stopped within a tactical approach conditional upon the events of the mission (grades and land forms, the supposed enemy position, risk-taking, orders, friendly position …) and it is subjected to a tactical move involving a progress by successive jumps.

The interest of separating the operating sequences of the vehicle is also stated clearly by the statistical characteristics of the power demand of the auxiliary consumers shown in Table III. Indeed, there is a noticeable difference in the consumption values depending on the vehicle operating mode. This remark suggests consequently considering three separate operating sequences: the first one when the vehicle is moving, the second and third one associated with the vehicle stop phases, with two further distinctions depending on whether the engine is switched off or not.

C. Cross-correlation analysis

The present section is now dealing with the analysis of the serial correlation for the onboard power demand of the auxiliary equipment. The aim is to detect any significant timedependence to be taken into account for further modelling purposes. The scatter graphs of Fig. 6, Fig 7 and Fig. 8 associated with the three operating sequences mentioned above are displayed both with the original measurements and the ranked values for the purpose of calculating Pearson and Spearman's rank correlation. A fairly strong relationship between the auxiliary powers for different time lags is pointed out at zero speed, particularly when the engine is switched off. Conversely, the operating sequence when the vehicle is running shows a far broader spread of power values at two distinct times. This view is confirmed by the Spearman-rank correlation coefficient assessments of Table IV which provides evidence there is a fairly strong link between past and current values of the power demand at zero-speed while crosscorrelations vanish outside of the stop sequences. V. CONCLUSION

The present study raises a specific problem of military vehicles within the use of auxiliary equipment and devices in terms of power consumption during their mission. At this point, the application of statistical methods to explore the possible correlations (inter-correlations and cross-correlations) between the auxiliary power and the parameters of the mission i.e. speed, acceleration, and road gradient is appropriate. The results lead then to conclude that the model of the electric consumption of the auxiliary equipment is to be reasonably separated from the model of the vehicle mission, because of a weak relationship pointed out by the Spearman's rank order correlation coefficient. Nevertheless, the data analysis shows the need for implementing three different auxiliary power models according to whether the vehicle is running or is stopped with the engine switched off or not. In this regard, the authors propose a methodology of analysis that can be easily reproduced for any other vehicles and that might bring in turn to different conclusions. Then, the interest of the approach is also to suggest directions for selecting suited models of the onboard power demand that would be able to take account of the cross-correlations put forward by the data analysis. Finally, it provides a basis for future research that would lead to a proper sizing of the electrical sources such as storage devices and generators.

Fig. 1 .

 1 Fig. 1. Recordings of one mission over a military training exercice: (a) Speed -(b) Acceleration -(c) Road gradient

Fig. 2 .Fig. 3 .

 23 Fig. 2. Vehicle speed and auxiliary power patterns over one full mission

Fig. 4 .

 4 Fig. 4. Scatterplot: speed versus auxiliary power for v > 0, (a) original measurements -(b) ranked values

Fig. 5 .

 5 Fig. 5. Auxiliary power demand at zero speed (a) engine off -(b) engine on

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Scatter graph: Paux(k) versus Paux(k-1) over stop sequences -engine 'off' (a) original measurements -(b) ranked values

TABLE I . GLOBAL MILITARY CHARACTERISTICS Parameter Unit Value

 I

		Distance	km	[24 -55]
		Duration	s	[15500 -47000]
	Speed	max	km/h [64 -74.4]
		average	km/h [6.4 -22.9]
	Acceleration	average	m/s 2	[-2.28 -1.384]
	Road gradient	%	[-20.76 -17.94]
		duration range	s	[2 -10500]
	Stops	number mean duration	--s	[92 -131] [313 -335]
	Equipment energy by mission	kwh	[23.44 -31.42]
	Energy mobility by mission	kwh	[110 -254]

TABLE II .

 II

CORRELATION COEFFICIENTS BETWEEN THE AUXILIARY POWER AND THE MISSION PARAMETERS.

Correlation coef. (%) Speed /Paux Acc. / Paux Slope / Paux Full mission (v ≥ 0)

		Pearson	34.02	-0.13	-1.68
		Spearman	44.13	4.19	-1.27
		Pearson	-12.83	-0.19	-4.64
	v > 0	Spearman	-7.31	2.34	-3.42

TABLE III .

 III

STATISTICAL CHARACTERISTICS OF THE AUXILIARY POWER DEPENDING ON THE VEHICLE SPEED Paux (

kW) full mission (v ≥ 0) vehicle running (v > 0) vehicle stopped (v = 0) regardless engine state engine 'off' engine 'on' mean value

		1.79	4.76	1.38	0.66	3.99
	standard deviation	2.21	2.70	1.79	0.40	1.92
	coef. of variation(%)	1.23	0.57	1.30	0.61	0.48