Alain Brenier

Andrzej Majchrowski

Edward Michalski

Light propagation properties of the Bi2ZnOB2O6 acentric biaxial crystal: angular orbital momentum from conical diffraction

may come

Light propagation properties of the Bi2ZnOB2O6 acentric biaxial crystal: angular orbital momentum from conical diffraction

Introduction

The acentric biaxial Bi2ZnOB2O6 (BZBO) oxy-borate crystal can be elaborated in bulk size by top-seeded Kyropoulos growth method [START_REF] Li | A Novel Nonlinear Optical Crystal Bi2ZnOB2O6[END_REF][START_REF] Li | Growth, Structure, and Optical Properties of a Congruent Melting Oxyborate, Bi2ZnOB2O6[END_REF]. It is non-hygroscopic and stable in a moist environment. Neodymium doping is also possible with high crystal quality [START_REF] Brenier | Chirality versus dichroism: competition and role in conical diffraction displayed from the Nd:Bi2ZnOB2O6 acentric biaxial laser crystal[END_REF]. Infra-red spectroscopy, transmittance and refractive indices measurements have been done. It was found that phase matching exists for second harmonic generation of the 1064 nm YAG:Nd laser radiation and its nonlinear optical coefficients have been determined by the Marker fringes method. So up to now this crystal received attention because it belongs to the wellknown family of the acentric borate crystals endowed with quadratic nonlinearities and promising optical properties: high transparency in the UV domain and a high optical damage threshold. They are ideal for responding to the need of high power UV and visible all-solid state light sources, up to the fifth harmonic generation (213 nm) of the 1064 nm YAG:Nd fundamental wavelength in some cases. Their optical properties depend on borate anionic groups in a comprehensive modelling [START_REF] Chen | Development of new nonlinear optical crystals in the borate series[END_REF][START_REF] Sasaki | Recent development of nonlinear optical borate crystals: key materials for generation of visible and UV light[END_REF][START_REF]Structure-Properties relationship in Non-Linear Optical Crystals I[END_REF].

On the other hand biaxial crystals, acentric or not, have an interesting property of linear optics: they can transfer an orbital angular momentum (OAM) to light propagating in directions close to the crystal optical axis. The OAM of an optical beam is the effect of the linear momentum acting off-axis with respect to its centre [START_REF] Padgett | Orbital angular momentum 25 years on[END_REF]. Helical phase-front beams, for example Laguerre-Gaussian modes, have such 𝑙 ℏ/photon OAM related to their exp (-𝑖𝑙𝜑) transverse phase variation. This momentum is due to the wave structure and can be in addition to the spin angular momentum (SAM) which is due to the polarization. Several devices can be used to generate beams with OAM. Let us cite the insertion in the beam path of a phase-plate with a thickness increasing with the azimuthal angle [START_REF] Beijersbergen | Helical-wavefront laser beams produced with a spiral phaseplate[END_REF] or a spatial light modulator based on liquid crystals whose refractive index is modified by an electric field [START_REF] Jesacher | Wavefront correction of spatial light modulators using an optical vortex image[END_REF]. Insertion of cylindrical lenses is another route [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF].

Going back to biaxial crystals, the propagation along the optical axis of a focused beam leads to its conical refraction (CR, instead of the usual double-refraction) and so its emergence as a hollow cylinder [START_REF] Landau | Electrodynamique des milieux continus[END_REF][START_REF] Born | Principles of optics: electromagnetic theory of propagation, interference and diffraction of light[END_REF][START_REF] Turpin | Conical refraction: fundamentals and applications[END_REF]. Nowadays applications of CR have been found, let us only cite the optical trapping of light [START_REF] O'dwyer | Conical diffraction of linearly polarized light controls the angular position of a microscopic object[END_REF]. The key-point for the present work is that the beam emerges also with modified OAM and SAM as it was calculated by Berry [START_REF] Berry | Orbital and spin angular momentum in conical diffraction[END_REF]. Theoretical predictions were nicely experimentally verified in the case of the centro-symmetric KGd(WO2)4 crystals [START_REF] O'dwyer | Generation of continuously tunable fractional optical angular momentum using internal conical diffraction[END_REF]. This work is devoted to the biaxial BZBO crystal. Because it is acentric, the light propagation is modified by the optical activity (this is a Pasteur medium) as we will show. We experimentally exhibit the CR of a laser beam in the near and far fields, and we study the intensity and the phase of the emergent waves. The connection with the OAM and SAM is obtained through a theoretical model of wave propagation including the birefringence, the optical activity and more generally the crystal bi-anisotropy.

Light propagation through a bi-anisotropic material

Far field pattern

Let us first find the plane wave eigen-modes propagating through a bi-anisotropic material.

Their complex electric field is written as 𝐄(𝐫)exp(𝐤 𝑡 ′ . 𝐫) with a similar writing for 𝐇, B and D fields and where 𝐤 𝑡 ′ has the meaning of the refracted (and transmitted) wave-vector inside the material. We can choose a transverse frame x1, x2, x3 whose third axis x3 is parallel to 𝐤 𝑡 ′ (𝑟, 𝜑), (𝑟, 𝜑) being the polar and azimuthal angles defined from the x', y', z' frame in Fig. 1 (z' is the optical axis of the biaxial crystal). More details can be found in ref. [START_REF] Brenier | Lasing with conical diffraction feature in the KGd(WO4)2:Nd biaxial crystal[END_REF]. Inserting these field expressions inside Maxwell-Ampère and Maxwell-Faraday equations, we find the following relations between the components:

[ 𝐸 1 𝐸 2 ] = 𝜔 k 𝑡 ′ [ 0 1 -1 0 ] [ 𝐵 1 𝐵 2 ] (1)
[ 𝐻 1 𝐻 2 ] = - 𝜔 k 𝑡 ′ [ 0 1 -1 0 ] [ 𝐷 1 𝐷 2 ] (2)
On the other hand, the constitutive relations of a bi-anisotropic material [18][19] in the DB representation are:

𝐄 = 𝛋. 𝐃 + 𝛘. 𝐁 (3)
𝐇 = 𝛄. 𝐃 + 𝛎. 𝐁 (4)
At this step the 𝛋, 𝛘 , 𝛄 and 𝛎 are 3-D tensors with components written in any frame. Then the components in Eq. ( 3) and ( 4) can be transformed towards the (x1, x2, x3) transverse frame (with an underlined notation) from their expression known in the dielectric xyz frame. For example the inverse permittivity tensor is expressed as 𝛋:

𝛋(𝑟, 𝜑) = 𝑆 2 (𝑟, 𝜑)(𝑆 1 𝛋𝑆 1 -1 )𝑆 2 -1 (𝑟, 𝜑) (5)
where 𝑆 1 is the 3-D rotation making the (x, y, z)(x', y', z') frame transfer and 𝑆 2 (𝑟, 𝜑)

makes the (x', y', z')( x1, x2, x3) transfer. The DB representation is very convenient because the D and B fields are transverse so D3=B3=0 and Eq. ( 3) and ( 4) simplify as:

[ E 1 E 2 ] = [κ] 2X2 [ D 1 D 2 ] + [χ] 2X2 [ B 1 B 2 ] (6)
[ H 1 H 2 ] = [𝛾] 2X2 [ D 1 D 2 ] + [ν] 2X2 [ B 1 B 2 ] (7)
where 2X2 means the 2-D restriction of the tensors to the 1 and 2 transverse coordinates.

We can determine the four transverse fields of each eigen-mode by combining the four (1),

(2), ( 6) and ( 7) equations. With Eq. ( 2) introduced inside Eq. ( 7) we obtain:

[ B 1 B 2 ] = [[ν] 2X2 ] -1 {-[γ] 2X2 + 𝜔 k 𝑡 ′ [ 0 -1 1 0 ]} [ D 1 D 2 ] (8)
Then with Eq. ( 1), ( 6) and ( 8) we find the 𝐃 eigen-modes as the non-vanishing solutions of:

{[κ] 2X2 -𝜇 0 [χ] 2X2 [𝛾] 2X2 + 𝜇 0 𝜔 k 𝑡 ′ {[χ] 2X2 [ 0 -1 1 0 ] -[ 0 -1 1 0 ] [γ] 2X2 } - 𝜇 0 ( 𝜔 k 𝑡 ′ ) 2 [ 1 0 0 1 ]} [ D 1 D 2 ] = [ 0 0 ] (9)
where we have restricted to the case useful in the present work of an isotropic and nonmagnetic permeability 𝛎 = 1/𝜇 0 .

Equating to zero the determinant leads to the generalisation of the well-known Fresnel equation which is a fourth degree equation including odd (1 and 3) degrees. Its two solutions in a given propagation direction are the eigen-values

𝜔 k 𝑡± ′ leading to the two [ D 1 D 2 ]
± eigen-modes. The eigen-modes of the other fields and specially the electric field are obtained from Eq. ( 1), ( 2) and ( 8).

The electric field 𝐄 of the input beam at the entrance face of the crystal is a priori known in the (x', y', z') frame. This field can be decomposed in plane waves from its 2D-Fourier transform 𝐄 ̂. Each plane wave propagates inside the crystal after refraction according to a calculation detailed in Ref. [START_REF] Brenier | Lasing with conical diffraction feature in the KGd(WO4)2:Nd biaxial crystal[END_REF] and summarized as:

-First step: we obtain the 𝐄 ̂components in the (x1, x2, x3) transverse frame -Second step: with the help of the eigen-modes we obtain the 𝐄 ̂tangential components at the output crystal face after a path L.

-Third step: we go back to the (x', y', z') frame which leads to the [𝐄 ̂(𝐿, 𝑘 𝑖𝑥 ′ , 𝑘 𝑖𝑦 ′ )] 𝑥′𝑦′ electric far field pattern.

Near field pattern after the material

After propagation in air over a distance d behind the crystal, all the plane waves constituting the far field interfere, in other words the electric near field is obtained as an inverse 2D-

Fourier transform [𝐄(𝑥 ′ , 𝑦 ′ , 𝑑)] 𝑥′𝑦′ = 𝑖𝐹𝑇 {exp(𝑖𝑘 𝑖𝑧 ′ 𝑑) [𝐄 ̂(𝐿, 𝑘 𝑖𝑥 ′ , 𝑘 𝑖𝑦 ′ )] 𝑥 ′ 𝑦 ′ }.
From discretisation the inverse Fourier transform is calculated with the inverse Fast Fourier Transform (iFFT) algorithm of the Mathlab package.

OAM and SAM after the material

In the paraxial approximation used here for calculating the light propagation in air after the conical diffraction inside the crystal, the beam OAM component along the propagation direction involves only the transverse electric field. The OAM/photon is obtained by the expression [START_REF] Berry | Orbital and spin angular momentum in conical diffraction[END_REF]:

𝐽 𝑜𝑟𝑏 = ℏIm ∬ 𝑑𝑹(𝐄 * .𝜕 𝜑 𝐄)
∬ 𝑑𝑹(𝐄 * .𝐄) [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] while the SAM/photon is given by: 𝐽 𝑠𝑝 = ℏIm ∬ 𝑑𝑹(𝒆 𝑧 .𝐄 * X𝐄) ∬ 𝑑𝑹(𝐄 * .𝐄) [START_REF] Landau | Electrodynamique des milieux continus[END_REF] These formulas show that the OAM depends on the angular variation 𝜑 of the electric field around the centre of the beam and that the SAM depends on the polarization.

Experimental set-up for CR and OAM visualization

BZBO single crystals were grown from stoichiometric melts by means of Kyropoulos technique on [001] oriented seeds. Detailed description of BZBO crystallization and morphology of as-grown single crystals can be found in our previous paper [START_REF] Brenier | Chirality versus dichroism: competition and role in conical diffraction displayed from the Nd:Bi2ZnOB2O6 acentric biaxial laser crystal[END_REF].

The beam from a 632.8 nm Helium-Neon laser, focused into a 16 µm waist spot behind the BZBO sample (1 cm) with a 5 cm focal length lens L1, is schematically represented in Fig. 1. The intensity of the near field (red circle) resulting of conical diffraction behind the crystal is imaged through the L2 lens (10 cm focal length) with a magnification 15 on the ocean Optics CCD screen, the magnification being calculated from the lens-camera distance. The far field is obtained in a separate measurement replacing L2 by a L3 lens (20 cm focal length and not shown in the figure), the CCD screen being located in the L3 focal plane. All the lenses are achromatic doublets from Thorlabs. In Fig. 1 we can see that a linear polarizer and an achromatic quarter-wave plate (from Thorlabs) can be introduced in the beam path before the sample in order to make the beam polarized linearly, circularly or elliptically. After the sample similar optical components are used to analyse the beam and decompose it in different polarization states.

In order to visualise the phase of the near field, a reference beam is taken from the input beam with a 50/50 beam splitter, is expanded by a telescope and redirected towards the CCD screen by a second 50/50 beam splitter. The whole set-up is similar to a Mach-Zehnder interferometer and the two beams go through the same output linear polarizer Pol2 allowing the interference. On the output beam-splitter face the two beams are a few mm separated (horizontally, vertically or in any other direction) in order to create a controlled wedge fringe pattern with the desired resolution.

Results and interpretation

Characterisation of the optical activity of the Bi2ZnOB2O6 acentric crystal

The BZBO crystal is orthorhombic with point group mm2 (C2v). The x, y, z axes of the dielectric frame coincide respectively with the b, c, a crystallographic axes. So the optical axis of the transparent crystal lies in the xz=ba plane with the Vz angle with the z axis close to 42.6°, depending slightly on the wavelength. The inverse permittivity tensor 𝛋 is diagonal in the xyz frame and its components are real (we limit this work to the domain of transparency of the material) and obtained through the measured three principal refractive indices [START_REF] Li | Growth, Structure, and Optical Properties of a Congruent Melting Oxyborate, Bi2ZnOB2O6[END_REF]: n1=2.0679, n2=2.105, n3=2.1541 .

The 𝛘 and 𝛄 tensors are constituted of a symmetric and an antisymmetric part. However in the present work our goal is to describe the optical activity as we show in the following of this sub-section. In this case only the antisymmetric part is necessary:

𝛘 = 𝝌 𝑎 (12) 𝛄 = -𝝌 𝑎 𝑇 ( 13
)
where the upper-script T means the matrix transpose. In this case the generalized Fresnel equation is biquadratic (no odd degrees). In the mm2 BZBO symmetry there are only two independent parameters, real in the domain of transparency:

where the tensor components are written in the xyz dielectric frame. The two 𝜒 𝑎𝑥𝑧 and 𝜒 𝑎𝑧𝑥 can be investigated as follow.

In this sub-section, removing the two quarter-wave plates in the experimental set-up we can launch a beam with a vertical polarization and detect it on the CCD with a horizontal analyser.

As a preliminary step the BZBO sample is replaced by a centro-symmetric crystal with the same orientation (z') of its optical axis, let us say KGd(WO2)4 (KGW). The result of the far field intensity is shown in Fig. 2 (a1) (experimental) and (a2) (theoretical). We can see that the centre is extinguished because it corresponds to propagation through an isotropic direction between two crossed polarizers.

To the contrary the similar experiment conducted with the BZBO sample leads to a brilliant centre in Fig. 2 (a3). The reason is that the initial polarization has been rotated by an angle measured to be 100° (it is easily measured by searching with the output polarizer the angle giving the light extinction). The specific rotatory power 𝜌 (rad/cm) is due to the fact that the two eigen-modes are left and right circularly polarized with a difference ∆𝑘 in their wavenumbers:

𝜌 = ∆𝑘/2 (15)
We measured 100° rotation for 0.6375 cm sample thickness, so =2.74 rad/cm. Because ∆𝑘 is related to the discriminant of Fresnel equation we find the relation:

χ 𝑎11 + χ 𝑎22 = 𝜆𝑐 𝜋𝑛 𝑦 2 𝜌 (16)
The two independent parameters in Eq. ( 14) cannot be determined separately from this experiment, only their sum is determined:

𝜒 𝑎𝑥𝑧 + 𝜒 𝑎𝑧𝑥 = - χ 𝑎11 +χ 𝑎22 cos(𝑉 𝑧 )sin (𝑉 𝑧 ) (17)
In the present work we will not conduct more refined experiments to determine separately these two parameters, but we will show that using as a try the same value for both of them half the value from Eq. ( 17), we obtain a satisfactory description of all the experimental figures. The theoretical Fig. 2 (a4) can be compared to Fig. 2

(a3).

The optical activity has also a noticeable influence in the phase of the wave despite that in this sub-section the input and output beams have no OAM. The intensity of the near field in the focal plane is represented in Fig. 2 (b1) (experimental) and 2 (b2) (theoretical) in the case of the centro-symmetric KGW. The intensity vanishes in the horizontal xz-plane which constitutes a border such that the phase of the wave above and below have an abrupt  rad difference as it is shown in Fig. 2 (c1) (experimental) and 2 (c2) theoretical) following the brilliant fringes. To the contrary in the BZBO case, the optical activity not only modifies the field intensity (Fig. 2 (b3) and (b4)) but also the phase: we can see following through the horizontal xz-plane the brilliant fringes in Fig. 2 (c3) and (c4) that they are no more head to foot with the dark fringes: the phase variation is smoothed.

Configuration generating an angular orbital momentum close to integer

In this sub-section, thanks to a QW plate and a vertical polarizer in the experimental set-up we launch a beam through the BZBO sample with a left-circular polarization with the goal to project the output beam on a right-circular state. For that purpose we detect the output beam on the CCD through a right-circular polarizer constituted with a second QW with its fast axis perpendicular to the one of the input QW plate, followed by a horizontal analyser and finally followed by a third QW plate. However the third QW plate is not actually installed because it will change neither the intensity nor the phase of the output wave and it should prevent the interference with the reference wave. In the focal plane we observe mainly two brilliant rings separated by a dark one (Fig. 3 (a1) experimental) whose diameter is 240 µm, that is to say much higher than the beam diameter (32 µm). The experimental diameter is very close to the calculated one (Fig. 3 (a2)). The interference pattern with the reference beam is shown in Fig. 3 (b1). We can see the modulation of the two brilliant rings by the alternation of dark and brilliant fringes. The keypoint is that there are 19 brilliant fringes in the upper external semi-circle and 20 in the lower one. Similarly there is one brilliant fringe more in the lower semi-circle of the internal ring.

The meaning of the surplus fringe if that the right circular output wave includes a 2 phase variation in addition to the phase of the reference wave and an inspection of the two optical paths of the waves (Fig. 1) allows to determine that this phase variation is incremental if we perform an anti-clock wise rotation around the centre of the circles. These experimental observations are in full agreement with the phase of the right circular output wave calculated by the model of section 2 and represented in Fig. 3 (b2). More, the OAM calculated with Eq. ( 10) is found to be 0.996 ℏ/photon and the SAM is -1 ℏ/photon. Let us add that the dark ring is a boundary on which the phase is not defined (Fig. 3 (b2)). It separates the field in two domains with a  rad phase difference visualised by head to foot brilliant and dark fringes in Fig. 3 (b1). Fig. 3 (b2) exhibits also an optical vortex in its centre where the phase is not defined.

In the plane about 2 cm further (the Raman spot plane) the left-circular polarized component of the output wave collapse into a brilliant centre (not shown here). But the output right circular component beam collapses into a narrow brilliant ring with a much weaker diameter:

31 µm (Fig. 3 (c1) experimental and 3 (c2) theoretical). A similar description than in the ring plane can be performed: interference pattern with the reference beam (Fig. 3 (d1)) and calculated phase (Fig. 3 (d2)) with an optical vortex in the centre.

Configuration generating a fractional angular orbital momentum

Generally speaking the CR of an elliptical polarized beam leads to an output beam, after projection on another elliptical state, with a fractional OAM. In this subsection we launch a left-circular polarization beam through the BZBO sample and we project the output beam on a linear polarized state with -45° angle from the horizontal axis. The intensity in the focal plane is represented in Fig. 4 (a1) (experimental) and 4 (a2) (theoretical). The light is localised on the upper semi-circle with the presence of the dark ring. The phase of the wave is studied from the interference pattern with the reference beam (Fig. 4 (b1)). Following a closed path inside the dark ring as the one visualised by the green ellipse, and counting the number of brilliant fringes, we find that the phase variation is 0 along this round trip. This is confirmed by following the same path on the theoretical Fig. 4 (b2). To the contrary, following a closed path outside the dark ring as the red circle in Fig. 4 (b1), we find that there is one surplus fringe in the lower semi-circle corresponding to 2phase variation on a round trip. This is confirmed by following a similar phase trip in the theoretical Fig. The OAM calculated with Eq. ( 10) is found to be 0.456 ℏ/photon (and the SAM is 0 ℏ/photon because the polarization is linear).

In the plane about 2 cm further (the Raman spot plane) the beam intensity distribution is modified and we can see a brilliant area surrounded by dark ones: Fig. 4 (c1) (experimental) and 4 (c2) (theoretical). The experimental interference pattern with the reference beam is visualized in Fig. 4 (d1). We have drawn three small red ellipses. Each of them surrounds a field of fringes such that there is a surplus brilliant one either in the left part of the ellipse either in the right part. In the first case (one ellipse) the meaning is 2 phase variation in the clock-wise round trip and in the second case (two ellipses) 2 phase variation in the anticlockwise round trip. We have drawn a small green ellipse surrounding an area with no surplus brilliant fringe, meaning no phase variation following a round trip. Three white arrows have been added to link the centres of the three red ellipses with the three optical vortices found in the theoretical description of the phase represented in Fig. 4 (d2). A fourth white arrow has been drawn to link the centre of the green circle to the area calculated to be regular. A full agreement is observed between the calculated phase variation and the interference pattern. A large green circle is also drawn in Fig. 4 (d1) with no surplus brilliant fringe following a round trip. This is explained because it surrounds one clockwise and one anticlock-wise 2 phase variation areas. A large red circle is drawn with one surplus brilliant fringe following a round trip. This is explained because it surrounds one clockwise and two anticlock-wise 2

Conclusion

We have grown an acentric biaxial BZBO single crystal by the Kyropoulos method under conditions of low temperature gradients. It was oriented and cut perpendicular to the optical axis located in the ab plane at 42.6° from the a-axis. We focused slightly behind the sample the Gaussian beam from a He-Ne laser at 632.8 nm and we made magnified images of the near field on a CCD camera. Linear polarizers and quarter-wave plates were used to launch and detect beams with various circular, elliptic or linear polarization states.

A reference beam was taken from the entrance and redirected on the CCD screen to obtain a controlled wedge fringe pattern visualising the phase of the wave. The specific rotatory power was measured to be 2.74 rad/cm. We present results of conoscopy patterns obtained theoretical.

Fig. 1

 1 Fig. 1 Experimental set-up for CR and OAM visualisation.

Fig. 2

 2 Fig. 2 See text in subsection 4.1. Column a: intensity of the far field; column b: intensity of the near field in the focal plane; column c: 1 and 3: interference with the reference beam, 2 and 4: phase in the focal plane (vertical scale in rad.). Lines 1 and 3: experimental, lines 2 and 4: theoretical.

Fig. 3

 3 Fig. 3 See text in subsection 4.2. Columns a and b: focal plane, columns c and d: Raman spot plane. Columns a and c: near field intensity; column b and d: 1: interference with the reference beam and phase (vertical scale in rad.). Line 1: experimental, lines 2: theoretical.

Fig. 4

 4 Fig. 4 See text in subsection 4.3. Columns a and b: focal plane, columns c and d: Raman spot plane. Columns a and c: near field intensity; column b and d: 1: interference with the reference beam and phase (vertical scale in rad.). Line 1: experimental, lines 2: theoretical.

  . Once more these phase variations can be found also following similar paths in the theoretical Fig.4(d2).

  launching through the crystal slab a left-circular polarized beam and projecting the output beam on a right-circular polarization state or on a linearly polarization state at -45° from the horizontal axis. All the experimental results (optical activity, intensity and phase of the output wave and vortices) are well described by a theoretical model based on the bianisotropy of the BZBO crystal. The OAM and SAM are calculated for each configuration from the distribution of the transverse electric field. 18. Jin Au Kong, Electromagnetic wave theory, EWMH publishing, Cambridge, Massachusetts, USA. 19. M. V. Berry, "The optical singularities of bianisotropic crystals", Proc. R. Soc. A, 461, 2071-2098 (2005).

Figure captions 1 . 3 . 4 .

 134 Figure captions1. Experimental set-up for CR and OAM visualisation. 2. column a: intensity of the far field; column b: intensity of the near field in the focal plane; column c: 1 and 3 interference with the reference beam, 2 and 4 phase in the focal plane (vertical scale in rad.). Lines 1 and 3: experimental, lines 2 and 4: theoretical. 3. See text in subsection 4.2. Columns a and b: focal plane, columns c and d: Raman spot plane. Columns a and c: near field intensity; column b and d: 1: interference with the reference beam and phase (vertical scale in rad.). Line 1: experimental, lines 2: theoretical. 4. See text in subsection 4.3. Columns a and b: focal plane, columns c and d: Raman spot plane. Columns a and c: near field intensity; column b and d: 1: interference with the reference beam and phase (vertical scale in rad.). Line 1: experimental, lines 2:


D Maamria

K Gillet

G Colin

Y Chamaillard

C Nouillant

Optimal predictive eco-driving cycles for conventional, electric and hybrid electric cars

In this paper, the computation of eco-driving cycles for electric, conventional and hybrid vehicles using receding horizon and optimal control is studied. The problem is formulated as consecutive-optimization problems aiming at minimizing the vehicle energy consumption under traffic and speed constraints. The impact of the look-ahead distance and the optimization frequency on the optimal speed computation is studied to find a trade-off between the optimality and the computation time of the algorithm. For the three architectures considered, simulation results show that in urban driving conditions, a look-ahead distance of 300m to 500m leads to a sub-optimality less than 1% in the energy consumption compared to the global solution. For highway driving conditions, a look-ahead distance of 1km to 1.5km leads to a sub-optimality less than 2% compared to the global solution.

I. INTRODUCTION

Nowadays, energy efficiency is increasingly becoming a major concern, and a subject of attention from major international organizations. In the automotive industry, the main research directions towards improving energy efficiency focus on decreasing carbon dioxide (CO 2 ) emissions. Engine technology and car performances have improved rapidly thanks to alternative propulsion technologies (hybrid and electric systems) and engine downsizing, while drivers have not adapted their behavior (driving style). Eco-driving can be defined as a multi-criteria optimization problem (fuel consumption, trip duration, drivability, etc.) of various tasks (navigation, guidance, stabilization) under safety constraints [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF]- [START_REF] Kim | Development of ecodriving guide system[END_REF]. In other words, the idea of eco-driving is to compute the vehicle speed and gear-box ratio set points that minimize the vehicle energy consumption. The results reported in the literature show the dependance of the energy saving via ecodriving on the nature of the trip. In urban driving conditions, the energy saving may reach 20% while for highway trips, the energy saving is about 5 to 10% [START_REF] Kim | Development of ecodriving guide system[END_REF]- [START_REF] Daun | Evaluation of driving behavior and the efficacy of a predictive eco-driving assistance system for heavy commercial vehicles in a driving simulator experiment[END_REF].

To help the drivers to reduce the energy consumption of their vehicles, a human machine interface (HMI) integrated in the dashboard proposes the optimal speed and gear-box ratio set points to the driver. The computation of these setpoints can be formulated as an optimization problem under constraints such as speed limits, traffic information, the trip duration and the final destination [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF], [START_REF] Ozatay | Analytical solution to the minimum energy consumption based velocity profile optimization problem with variable road grade ?[END_REF]- [START_REF] Mensing | Ecodriving: an economic or ecologic driving style?[END_REF]. Depending on the horizon of the available information, two scenarios are distinguished: Off-line scenario and On-line scenario.

In the off-line scenario, the constraints (speed limits and traffic information) are fully known in advance. This question was addressed for conventional vehicles in [START_REF] Mensing | Optimal energy utilization in conventional, electric and hybrid vehicles and its application to eco-driving[END_REF], [START_REF] Maamria | Optimal eco-driving for conventional vehicles: simulation and experiment[END_REF], [START_REF] Sciarretta | Optimal ecodriving control: Energy-efficient driving of road vehicles as an optimal control problem[END_REF], for electric cars in [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF], [START_REF] Dib | Optimal energy management for an electric vehicle in eco-driving applications[END_REF], [START_REF] Miyatake | Theoretical study on ecodriving technique for an electric vehicle considering traffic signals[END_REF], [START_REF] Maamria | On the use of dynamic programming in eco-driving cycle computation for electric vehicles[END_REF] and for hybrid electric cars in [START_REF] Hellstrom | Design of an efficient algorithm for fuel-optimal look-ahead control[END_REF], [START_REF] Sciarretta | Optimal ecodriving control: Energy-efficient driving of road vehicles as an optimal control problem[END_REF], [START_REF] Bouvier | Determination and comparison of optimal eco-driving cycles for hybrid electric vehicles[END_REF]- [START_REF] Van Keulen | Velocity trajectory optimization in hybrid electric trucks[END_REF]. In this case, an upper bound on energy saving that can be reached through eco-driving, can be computed and used as a reference. In addition, some heuristic rules to reduce the fuel consumption may be extracted from the optimal solutions: anticipate the traffic flow, drive at a low engine speed. The ability to anticipate traffic flow depends on many factors such as the road layout, visibility and it is therefore unreasonable to expect it consistently of a driver. It is necessary to use algorithms to gather data from available sources and help the driver to anticipate.

In the on-line scenario, the future driving conditions are partially known for a given time or distance horizon usually called e-horizon. The information can be retrieved and updated in real time from many sources such as onboard sensors, ADAS (Advanced Driver-Assistance Systems) and navigation systems (Global Positioning System GPS unit). Based on the available information, the optimizer has to find the set points of the vehicle speed and the gear-box ratio in real time to minimize the energy consumption of the vehicle. This question has mainly been studied for heavy duty vehicles (trucks) due to the cost of fuel consumption for haulage companies [START_REF] Hellstrom | Look-ahead control for heavy trucks to minimize trip time and fuel consumption[END_REF]- [START_REF] Neiss | Predictive speed control for a motor vehicle[END_REF] where the roadway grade is among the biggest contributors to high fuel consumption. The objective was to find the speed and gear-box ratio set-points that minimize a trade-off between the fuel consumption and the trip duration on a predefined speed range. The studies in [START_REF] Schwickart | An efficient nonlinear model-predictive eco-cruise control for electric vehicles[END_REF], [START_REF] Kamal | Eco-driving using real-time optimization[END_REF] deal with an eco-driving system using Model Predictive Control (MPC) within a given prediction horizon for passenger cars in an urban road network with complex traffic flows. The study in [START_REF] Levermore | Real-time velocity optimization to minimize energy use in passenger vehicles[END_REF] highlights the potential of realtime vehicle speed optimization as a mechanism for energy saving. A driver feedback system, designed using Dynamic Programming (DP), was demonstrated in a test vehicle that uses the optimal velocity profile for the upcoming road. The DP algorithm was implemented with a distance step interval of 50m and 30 steps in the horizon, giving a total of 1.5km. The distance step interval of 50m from our viewpoint is high (the time needed to travel 50m when the vehicle speed is 100km/h is 1.8s) and the choice of a 1.5km horizon was chosen arbitrarily.

This paper follows the same approach as the on-line scenario and pushes it further for passenger (conventional, electric and parallel hybrid electric) cars. Depending on the length of the e-horizon, the on-line scenario will lead to a sub-optimality in the energy saving with respect to the global optimum (full knowledge of the future driving conditions). By using MPC and DP, the objective is to select the value of the look-ahead horizon in order to minimize the induced suboptimality in the energy saving while keeping the algorithm computation time reasonable. This choice is also interesting to size/choose sensors that predict future driving conditions. For this, various values of the e-horizon are considered and the induced sub-optimality compared to the global optimum is quantified. The choice of the look-ahead distance will be realized based on the optimality/computation time balance.

The paper is organized as follows. In Section II, the vehicle model and the computation of eco-driving cycles for conventional and electric cars are detailed. Section III presents the proposed predictive scheme. Numerical and simulation results are discussed in Section IV. An extension to hybrid electric cars is discussed in Section V. Some conclusions about the algorithm tuning are drawn in light of the presented results.

II. PROBLEM FORMULATION

A. Vehicle Modeling [START_REF] Guzzella | Vehicle propulsion systems[END_REF] 1) Motion equations: The vehicle is modeled in a vertical plane. According to Newton's law of motion, the vehicle speed v satisfies the following differential equation

m • dv(t) dt = F t (t) -F r (t), (1)
where m is the total vehicle mass including rotating parts, F t is the traction force provided by the prime mover and F r is the sum of resistance forces including rolling resistance force and aerodynamic drag force. Its expression is F r (t 2 , where c i , i = {0, 1, 2} are the coefficients of the road load equation.

) = c 0 + c 1 • v(t) + c 2 • v(t)
2) Prime mover: In the case of a conventional vehicle, the prime mover is an Internal Combustion Engine (ICE). In this study, the ICE is a Diesel engine. Its measured fuel consumption ṁ f (g/s) is computed through a look-up table as a function of the engine torque T e and the engine speed ω e : ṁ f (ω e , T e ) (see Figure 1) 1 .

The choice of a Diesel engine is not restrictive as the proposed method can be easily extended to gasoline engines.

In the case of an electric car, the prime mover is an electric machine. It is modeled by a quasi-static map describing the electric power (Figure 2) as a function of the electric machine torque T e and speed ω e : P m (ω e , T e ). This measured map usually includes the losses in the electric machine and the power electronic devices.

3) Transmission: The torque T e is related to the driver's torque demand at the wheel T wh by

T wh (t) = r tire • F t (t) = η gb • R gb (t) • R t • T e (t), (2)
1 The relation between ṁ f (g/s) and SFC (g/kWh) is given by where R gb is the gear-box ratio, η gb is the gear-box efficiency (assumed to be constant) and R t is the differential ratio. The torque T wh can be positive (traction) or negative (braking).

SFC = ṁ f 3600 × ω e × T e
Similarly, the rotational speed ω e is related to v by

ω e (t) = R gb (t) • R t • v(t) r tire . (3)
4) Battery model: The battery is represented by an equivalent circuit model comprising a voltage source U ocv in series with an electric resistance R b , both of which vary with ξ , the battery state of charge (SOC) [START_REF] Guzzella | Vehicle propulsion systems[END_REF], [START_REF] Badin | Hybrid Vehicles[END_REF]. The expression of the battery current I b is given by [START_REF] Guzzella | Vehicle propulsion systems[END_REF]

I b (ξ , v, T e ) = 1 2R b (ξ ) U ocv (ξ ) -U 2 ocv (ξ ) -4R b (ξ ) • P m (v, T e ) .
Based on the study [START_REF] Maamria | On the use of dynamic programming in eco-driving cycle computation for electric vehicles[END_REF], constant values of R b and U ocv are used: neglecting the dynamics of ξ in the eco-driving cycle computation leads to a negligible sub-optimality. The electrochemical battery power is defined by

P ech = I b •U ocv .
The dynamics of ξ is given by

dξ (t) dt = - I b (t) Q 0 ,
where Q 0 is the nominal battery capacity. The model parameters are summarized in Table I. The numerical values are omitted for confidentiality reasons.

B. OCP formulation

The eco-driving approach consists in finding the best vehicle speed profile that minimizes the vehicle energy consumption for a fixed road. The vehicle starts from a point A at a velocity v 0 (≥ 0) and must reach a destination point B at time t f , with a velocity v 1 (≥ 0). Some constraints have to be considered: acceleration and speed limits, and traffic constraints [START_REF] Dib | Optimal energy management for an electric vehicle in eco-driving applications[END_REF], [START_REF] Mensing | Ecodriving: an economic or ecologic driving style?[END_REF]. This question can be formulated as an Optimal Control Problem (OCP) [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF], [START_REF] Mensing | Ecodriving: an economic or ecologic driving style?[END_REF].

The cost function to be minimized is

J = t f 0 L(ω e (t), T e (t))dt,
where t f is the trip duration. In the case of a conventional vehicle, L is the fuel consumption ṁ f while in the case of an electric car, L is the electrochemical battery power P ech .

The control variable to be optimized u is composed of two components: the torque T e and the gear-box ratio R gb

u(t) = [T e (t), R gb (t)].
For the electric car, the considered gear-box has one ratio, only the electric machine torque is optimized. This optimization is carried out under the following stable dynamics

dv(t) dt = f (v(t), u(t)), v(0) = v 0 , (4)
dx(t) dt = v(t), x(0) = 0, ( 5
)
where x is the vehicle position and the non-linear function f is calculated by combining (1, 2)

f = 1 m (-c 0 -c 1 • v -c 2 • v 2 + η gb r tire • R gb • R t • T e ).
Since the speed, the torque T e and the gear-box ratio are limited, the final position and speed are fixed, the optimization is performed under the following mixed state and input constraints

v ∈ [0, v max (x)], (6) f (v, u) ∈ [a min , a max ], (7)
T e ∈ [T min (ω e ), T max (ω e )], (8
)

ω e ∈ [ω min , ω max ], (9) x
(t f ) = D, (10) v(t f ) = v 1 , (11)
where v max (x) is the speed limits at the vehicle position x, D is the total traveled distance, v 1 is the desired final speed, T min and T max are the minimum and the maximum prime mover torque given by look-up tables as a function of the rotation speed ω e (see Figures 1 and2). Equation ( 7) limits the vehicle acceleration between its maximum and minimum values. The acceleration is an algebraic function of the vehicle speed and the control variables. This function can be evaluated for all the possible choices of the vehicle speed and the control variables. The values not satisfying the constraints are excluded. The constraints on the prime mover torque in [START_REF] Ozatay | Analytical solution to the minimum energy consumption based velocity profile optimization problem with variable road grade ?[END_REF] and the rotational speed in ( 9) are mixed input-state constraints, they depend on the vehicle speed v and the gear-box ratio R gb .

In this paper, the torque converter dynamics is not considered in the OCP formulation because it increases the complexity of the optimization problem (thus the computation time) and the experimental and simulation results are very close [START_REF] Maamria | Optimal eco-driving for conventional vehicles: simulation and experiment[END_REF]. To summarize, the OCP is

(OCP) : min u t f 0 L(v, u)dt (12)
under the dynamics (4, 5), the state and input constraints [START_REF] Maamria | Optimal eco-driving for conventional vehicles: simulation and experiment[END_REF][START_REF] Daun | Evaluation of driving behavior and the efficacy of a predictive eco-driving assistance system for heavy commercial vehicles in a driving simulator experiment[END_REF][START_REF] Ozatay | Analytical solution to the minimum energy consumption based velocity profile optimization problem with variable road grade ?[END_REF][START_REF] Hellstrom | Design of an efficient algorithm for fuel-optimal look-ahead control[END_REF], and the final constraints [START_REF] Mensing | Ecodriving: an economic or ecologic driving style?[END_REF][START_REF] Sciarretta | Optimal ecodriving control: Energy-efficient driving of road vehicles as an optimal control problem[END_REF]. Thus, the optimized vehicle speed trajectory can be computed using (4) once the control variable u is optimized and the initial condition of the vehicle speed is known. 1) Speed limits: To specify the speed limits v max , a certain (fixed) margin e l (> 0) on the initial driving cycle speed is considered

v max (x) = v re f (x) + e l , v re f (x) > 0, 0, v re f (x) = 0,
where v re f (x) is the speed value of the initial driving cycle at the position x which is used as a reference of comparison to point out the benefit of eco-driving in terms of fuel consumption reduction. Note that the eco-driving cycle must have the same number of stops as the initial driving cycle [START_REF] Mensing | Ecodriving: an economic or ecologic driving style?[END_REF].

Other methods to compute speed limits can be considered [START_REF] Mensing | Optimal energy utilization in conventional, electric and hybrid vehicles and its application to eco-driving[END_REF], [START_REF] Bouvier | Determination and comparison of optimal eco-driving cycles for hybrid electric vehicles[END_REF]. For example, the speed limits can be given by the traffic signs or provided by a cloud server depending on the vehicle position. The choice of method used to compute the speed limits does not impact the solving method.

2) Solving method: The (OCP) defined in ( 12) can be solved using several methods: Pontryagin Minimum Principal (PMP) or Dynamic Programming (DP). Because of the non-linear nature of the optimization problem and the mixed state-input constraints, methods based on PMP are more complicated to implement [START_REF] Malisani | Dynamic control of energy in buildings using constrained optimal control by interior penalty[END_REF]. Thus, the method proposed here is based on Dynamic Programming [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF], [START_REF] Bryson | Applied optimal control[END_REF].

To reduce the computation time, the chosen approach transforms a time-based OCP into a distance-based OCP while introducing an additional term in the cost function as suggested in [START_REF] Monastyrsky | Rapid computation of optimal control for vehicles[END_REF], [START_REF] Mensing | Optimal energy utilization in conventional, electric and hybrid vehicles and its application to eco-driving[END_REF]. If the position space is discretized in N fixed steps of length ∆x, the time step ∆t(k), k = 1 : N is variable and is calculated from the vehicle speed v(k) and the acceleration a(k) by solving the following equation [START_REF] Mensing | Optimal energy utilization in conventional, electric and hybrid vehicles and its application to eco-driving[END_REF]

∆x = 1 2 a(k) • ∆t(k) 2 + v(k) • ∆t(k).
The final constraint on the vehicle position ( 10) is fulfilled by construction: the total traveled distance is D = N • ∆x. The trip duration is given by t f = ∑ N i=1 ∆t(i). To control t f , an additional tunable term β • ∆t(k) is added to the cost function [START_REF] Monastyrsky | Rapid computation of optimal control for vehicles[END_REF] as follows

Js (u) = k=N ∑ k=1 [L(v(k), u(k)) + β ] ∆t(k),
where the constant tunable parameter β penalizes the final time to obtain almost the same time duration as the initial driving cycle. The tuning of this parameter is discussed in the next paragraph.

3) Estimation of β : The parameter β is introduced to penalize the trip duration t f . This parameter is usually tuned iteratively [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF], [START_REF] Mensing | Optimal energy utilization in conventional, electric and hybrid vehicles and its application to eco-driving[END_REF]. The objective here is to estimate heuristically the value of β from the initial trip characteristics while keeping the error on the final time acceptable.

For this purpose, various driving cycles are used. The relation between β (searched iteratively) and the mean speed value v of the driving cycles (which is global information of the trip) is studied. The results in the case of the electric car for two values of e l are summarized in Figure 3. For a fixed value of e l , the relation between β and v is monotone: β increases when v increases. This heuristic relation can be approximated by a quadratic form (dashed curves in Figure 3) as follows

β = α 2 • v2 + α 1 • v + α 0 , (13)
where α i , i = {0, 1, 2} have to be identified. These parameters depend on the vehicle characteristics such as the weight and the coefficients of the road load equation. In addition, for the same value of v, β increases when e l decreases.

To validate the proposed relation to estimate β , the error on the final time for various driving cycles is given in Table II. The Artemis Highway [START_REF] Michel | The ARTEMIS European driving cycles for measuring car pollutant emissions[END_REF] and WLTC (the Worldwide harmonized Light vehicles Test Cycle) are two normalized driving cycles. The driving cycles 1 and 2 are real-driving cycles provided by the French car manufacturer PSA. The error on the final time is less than 3.9%, which is considered acceptable as the value of β is only estimated based on the mean speed value of the trip. The same approach is used for the conventional vehicle and the results are similar: the error on the final time is less than 4%.

The heuristic relation between β and the trip mean speed value simplifies the tuning procedure of β by providing a first good guess of its value and thus restricting the research field. This relation will be more useful for an embedded version of the proposed algorithm where the computation time and the number of iterations to find β have to be reduced.

III. PREDICTIVE ECO-DRIVING CYCLE COMPUTATION

It is well-known that control strategies obtained by optimal control use the full knowledge of the future (speed limits and traffic information in this study). However, these future conditions are uncertain and not totally known in advance in real situations. They may change in time and they depend on many factors such as local traffic, infrastructure status, non-vehicle actors and weather conditions. To deal with this situation, Model Predictive Control (MPC) techniques can be used [START_REF] Back | Predictive control of drivetrains[END_REF]- [START_REF] Debert | Predictive energy management for hybrid electric vehicles -prediction horizon and battery capacity sensitivity[END_REF].

A. Principle

The idea of the suggested method is the following: the future driving conditions are only known for a look-ahead distance D 0 ≤ D and the speed trajectory will be recalculated every traveled distance D f called optimization frequency where D f < D 0 . The total traveled distance D is divided into n intervals D i , i = 0 : n -1 where n is the ratio between D and D f . The OCP in (12) will be solved n times

(OCP i ) : min u t i+1 t i L(v, u)dt, (14)
under the dynamics (4, 5), the state and input constraints [START_REF] Maamria | Optimal eco-driving for conventional vehicles: simulation and experiment[END_REF][START_REF] Daun | Evaluation of driving behavior and the efficacy of a predictive eco-driving assistance system for heavy commercial vehicles in a driving simulator experiment[END_REF][START_REF] Ozatay | Analytical solution to the minimum energy consumption based velocity profile optimization problem with variable road grade ?[END_REF][START_REF] Hellstrom | Design of an efficient algorithm for fuel-optimal look-ahead control[END_REF]. The time t i corresponds to the moment when the vehicle reaches the position (i • D f ) while the time t i+1 corresponds to the moment for the position (i

• D f + D 0 ).
Thus, the constraints on the vehicle position are

x(t i ) = i • D f , x(t i+1 ) = i • D f + D 0 .
For i = n -1, the vehicle position will be (n -1)

• D f + D 0 .
If this position is greater than D, the look-ahead horizon will be only

D -(n -1) • D f ≤ D 0 .
As the future speed limits are only known for D 0 , vehicle speed is free at the end of each segment (i

• D f + D 0 ) except for the last segment v(t n ) = v 1 ,
where v 1 is the desired final speed at the end of the trip. Usually, this value is null, as it corresponds to a vehicle stop. This choice simplifies the algorithm as it is difficult to compute the appropriate speed at the position (i

• D f + D 0 ).
The optimization scheme is illustrated in Figure 4. The optimization problem for a total traveled distance of D is divided into n sub-problems solved successively. At each traveled distance D f , the speed limits and the traffic information are updated. A new vehicle speed trajectory is calculated on a distance horizon D 0 by solving the OCP (14) using the method described in Section (II-B.2). The speed trajectory obtained has to be followed by the driver. To avoid the effects due to human interaction, the vehicle is assumed to follow the speed set points precisely (in simulation).

B. Choice of D f

The optimality of the algorithm described in the previous section depends on the choice of D 0 and mainly on D f . The choice of D 0 depends on the information that can be gathered from the available sources (GPS and vehicle sensors). In this paper, for various values of D 0 , the choice of D f is optimized to find the best trade-off between the optimality and the implementation constraints of the algorithm.

• The value of D f must be greater than a lower bound D f ,min defined by the traveled distance needed to solve the optimization problem (14) (real-time computation constraint). The value of this distance depends on the value of D 0 , the equipment used (hardware) for computing and the available computation time: D f ,min increases when D 0 increases. D f ,min depends also on the vehicle speed: a low (respectively high) computation time for low (respectively high) mean speed value driving cycles. Thus, D f ,min will be chosen for highway driving cycles. In this study, D f ,min is 100m.

• An upper bound D f ,max on D f is computed based on the optimality. The methodology is the following: for a given value of D 0 and for a fixed trip duration t f , the value of D f is decreased from D 0 until D f ,min . The suboptimality is computed. Figure 5 shows the extra-fuel consumption and the variation of β for various values of D f when D 0 = 2km (β is tuned to get a constant trip duration): -Optimality is conserved for D f ≤ D f ,max : the solution is still optimal with a constant value of β . -Sub-optimality increases for D f > D f ,max : the suboptimality solution increases when D f tends to D 0 . In a distance window of length D 0 , the traffic constraints and speed limits are checked periodically: if there is a modification in the traffic constraints, the optimization is relaunched. By this choice, the use of the Electronic Control Unit (ECU) will be reduced as the optimization is performed only when it is necessary. Obviously, the optimization can be restarted every D f ,min , but the ECU will be constantly used. This behaviour increases the electric energy consumed by the ECU and thus increases the energy consumption.

0.5 1

IV. NUMERICAL ILLUSTRATION

To illustrate the proposed algorithm, five normalized driving cycles are used: ECE-15 (Urban Driving Cycle UDC), EUDC (Extra-urban driving cycle), WLTC (Worldwide harmonized Light vehicles Test Cycle), the Artemis Urban and Artemis Rural driving cycle [START_REF] Michel | The ARTEMIS European driving cycles for measuring car pollutant emissions[END_REF]. The duration t f without stop phases, the total traveled distance D, the final SOC (for the electric car) and the fuel consumption (for the conventional case) for each initial driving cycle (vehicle speed = v re f ) are given in Table III. The considered gearbox has 6 ratios for the conventional vehicle. The value of e l is arbitrarily fixed to 2km/h.

A standard computer equipped with an Intel Core i5-4590, CPU 3.30GHz with 8GB of RAM is used. Six scenarios are tested with the values of the look-ahead distance D 0 and the value of D f given in Table IV.

Scenario I II III IV V V I D 0 [km] D 3 2 1.5 1 0.5 D f [km] D
1.9 0.9 0.34 0.26 0.14

The values of D f were identified from a driving cycle of 6km with only one speed limit at 90km/h: an acceleration from 0km/h, a stabilized speed around 80km/h and a deceleration to 0km/h. The identification was done for the electric and conventional vehicles. The values of D f for the conventional vehicle were lower than the values of the electric case. The minimum value between the two cases is taken. Scenario I is considered as a reference of comparison (global optimum) in this study as the look-ahead distance D 0 equals the total traveled distance: full knowledge of the future driving constraints.

For the DP, a distance step of ∆x = 10m is used in the case of the ECE and the A. Urban and ∆x = 20m for the other driving cycles. These values of ∆x were chosen so as to have a time step approximately around 1s [START_REF] Maamria | Which methodology is more appropriate to solve eco-driving optimal control problem for conventional vehicles?[END_REF]. Two cases are studied: electric and conventional.

A. Case of an electric car

The cost function L is the electrochemical battery power P ech . The control variable is the electric machine torque T e . The following DP mesh parameters are used: ∆v = 0.02m/s for the vehicle speed and ∆T e = 2N.m For the electric machine torque. A sensitivity study using many driving cycles has been conducted to choose the best values of the quantization (∆v and ∆T e ). The criterion was to find a tradeoff between the solution optimality and its computation time. The initial value of the SOC is 90%.

For scenario I, the results in terms of energy consumption, energy saving through eco-driving with respect to the initial driving cycles and computation time t co are given in Table V. The extra-energy consumption and the error on the final time t f with respect to scenario I are given in Figure 6. For each cycle, the results are obtained using the same value of β calculated by [START_REF] Maamria | On the use of dynamic programming in eco-driving cycle computation for electric vehicles[END_REF]. For each scenario from Table IV, a positive value of the error on t f means that the trip duration obtained is greater than its value for scenario I. The consequence of that is a decrease in the energy consumption (as can be seen from the negative values of the extra-energy consumption in Figure 6). These results show that the trip duration t f increases slowly with a maximum error of 3s when D 0 decreases until 1km except for the Artemis urban case where the error is about 8s. When D 0 = 0.5km, the error on t f increases significantly: 3% of the total time in the WLTC case. At the same time, the energy consumption decreases.

To make a fair comparison between the various scenarios, the corrected extra-energy consumption with respect to scenario I is given in Figure 7. The correction of the energy consumption is calculated to take into account the difference in the trip duration as follows:

corrected energy consumption = energy consumption-β •∆t f ,
where ∆t f is the difference between the obtained final time and its target value. From Figure 7, the induced sub-optimality increases when the look ahead distance D 0 decreases, and its maximum value is less than 1.2%. This sub-optimality reduces the energy saving of eco-driving to 11.2% for the EUDC (see Table V). Two cases are distinguished:

I
• Driving cycles with a low mean speed value (ECE-15 and A. Urban): the sub-optimality is lower than 0.6%. • Driving cycles with a high mean speed value (A. Rural, WLTC and EUDC): scenario V , where sub-optimality is lower than 0.4%, is considered as a good choice. The vehicle speed trajectories for the EUDC with a constant value of β are given in Figure 8. For D 0 ≥ 1km, the speed trajectories are identical. When D 0 < 1km, drops appear in the vehicle speed trajectory because the vehicle speed is free at the end of each iteration as explained in Section III-A. That is why in some cases, the value of β needs to be adjusted. For this kind of cycle with high vehicle speed, the problem of speed drops appears when D 0 ≤ 1km even if D f = ∆x (not possible because of the implementation constraint on the choice of D f ). The mean computation time for each optimization and the number of optimizations for the various scenarios are given in Table VI. These values represent the mean time needed to get the new trajectory for a distance window of length D 0 . The computation time is approximately constant for each scenario (in the same column): the computation time per iteration for scenarios V (for the EUDC, the A. Rural and the WLTC) and V I (for the ECE-15 and the A. Urban) is 0.5s in the worst case. On the other hand, the overall computation time for the various scenarios increases (in the same row): to get the overall computation time, multiply the number of iterations by the mean time.

B. Case of a Conventional vehicle

The cost function L to be minimized is the fuel consumption ṁ f . The control variables are the engine torque T e and the gear-box ratio R gb . The following DP mesh parameters are chosen based on the results presented in [START_REF] Maamria | Which methodology is more appropriate to solve eco-driving optimal control problem for conventional vehicles?[END_REF]: ∆v = 0.1m/s for the vehicle speed and ∆T e = 2N.m for the engine torque2 . For scenario I, the results in terms of fuel consumption, fuel saving through eco-driving with respect to the initial driving cycles and computation time of the DP t co are given in Table VII. As for the electric case, to make a fair comparison between the various scenarios, the corrected extra-fuel consumption with respect to scenario I is given in Figure 9. The correction is introduced to make up the difference in the trip duration:

corrected fuel consumption = fuel consumption -β • ∆t f ,
where ∆t f is the difference between the obtained final time and its target value. The extra-fuel consumption induced by the predictive scheme reduces the fuel saving of eco-driving in Table VII: for scenario V I, the fuel saving with respect to the initial driving cycle is reduced to 14.5% in the case of the EUDC. The extra-fuel consumption increases when the look ahead distance D 0 decreases. In the case of the ECE-15 and the Artemis Urban, the maximum sub-optimality is less than 0.1% for all the scenarios considered. The sub-optimality reaches 0.4% for D 0 = 0.3km and D f = 0.1km for these two driving cycles. In the case of the Artemis Rural, the WLTC and the EUDC, the sub-optimality is less than 1% for scenario V .

The vehicle speed trajectories versus distance for the Artemis Urban are given in Figure 10. The trajectories are identical. The number of optimizations n and the mean computation time for each optimization for the various scenarios are given in Table VIII. The computation time for scenarios V (the EUDC, the A. Rural and the WLTC cases) and V I: the ECE-15 and the A. Urban cases) is less than 1.2s in the worst case (the time needed to travel 34m if the vehicle speed is 100km/h). These results for the electric and the conventional cases that it possible compute cycles in a reasonable time (1.2s in the worst case) with a partial knowledge of the future driving constraints and a quasioptimal energy consumption. These promising results are a first step toward the implementation of the eco-driving algorithm in an Engine Control Unit (ECU). A possible extension to hybrid cars is discussed in the next section.

V. EXTENSION TO THE HYBRID ELECTRIC VEHICLES A. Architecture description

The system considered here is a dual shaft parallel mild hybrid with an electric machine of 12kW connected to the Diesel engine described in Section II-A.2 by a belt (Figure 11). The gearbox is between the power-train and the wheel. Due to the choice of architecture:

• for a given gear-box ratio, the engine's speed and the electric machine's speed are dependent, • during the ZEV mode, the engine injection is cut off and the electric machine keeps the engine rotating.

For more details about this architecture, see [START_REF] Michel | Optimizing fuel consumption and pollutant emissions of gasoline-hev with catalytic converter[END_REF], [START_REF] Simon | Gasolinehev equivalent consumption and pollutant minimization strategy[END_REF]. For the transmission, the engine torque T eng and the electric machine torque T el are related to the torque required at the wheel T wh by

T wh (t) = η gb • R gb (t) • R t • [T eng (t) + R el • T el (t)] , (15)
where R el is the constant motor-to-wheel transmission ratio.

Similarly, the engine and the electric machine rotating speeds are linked by ω el (t) = R el • ω eng (t).

B. Eco-driving cycle computation: Off-line solution

The simplified approach suggested in [START_REF] Ngo | An optimal control-based algorithm for hybrid electric vehicle using preview route information[END_REF] for HEV is used as its induced sub-optimality with respect to the optimal solution is small (< 4%) while its computation time is reasonable. The method involves two steps:

1) Step1: An eco-driving cycle is calculated assuming that the vehicle is only propelled by the ICE (the vehicle parameters such as the weight and the road load coefficients do not change). The associated OCP is given in equation ( 12). 2) Step2: An energy management strategy (torque split) is calculated using dynamic programming. The vehicle is assumed to follow the eco-driving cycle obtained in Step1. The duration of the trip was fixed in the previous step. The control variable is the engine torque T eng .

As the torque T wh is imposed (the vehicle follows a pre-defined driving cycle with a predefined gear-box ratios), the electric machine torque T el is calculated from the torque balance in [START_REF] Kim | Two-stage optimal control of a parallel hybrid vehicle with traffic preview[END_REF]. A final constraint on the SOC is imposed. The formulation of this optimization problem was addressed in [START_REF] Guzzella | Vehicle propulsion systems[END_REF], [37] A similar approach was used in [START_REF] Sciarretta | Optimal ecodriving control: Energy-efficient driving of road vehicles as an optimal control problem[END_REF], [START_REF] Ngo | An optimal control-based algorithm for hybrid electric vehicle using preview route information[END_REF] where the second step was performed using the Pontryagin Minimum Principle (PMP) [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF].

As the two steps are decoupled, other methods to design the energy management strategy can be used: for example, the well-known ECMS (Equivalent Consumption Minimization Strategy) [START_REF] Paganelli | Equivalent consumption minimization strategy for parallel hybrid powertrains[END_REF]- [START_REF] Onori | and predictive control, optimization, identification, modeling for control, engine control, and energy management for hybrid vehicles[END_REF] or any other heuristic strategy.

C. Predictive Eco-driving cycle computation

In a first step, the method described in Section III-A to compute an eco-driving cycle for a conventional vehicle is used. Next, an energy management strategy is designed. For this, the driving cycle to be followed is only known for a look-ahead distance D 0 . The optimization will be done every traveled distance D f . The values of D 0 and D f given in Table IV are used.

As is commonly done for an HEV, a final constraint on the SOC is imposed in order to make a fair comparison between optimization strategies. A sustaining final constraint on the SOC will be imposed at the end of each segment:

ξ (t i+1 ) = ξ (t i ) = ξ re f , (16)
where ξ re f is the desired final SOC. Only the trajectory of ξ on a distance window of length D f is considered: the SOC does not come back to its target final value at the end of each traveled distance D f .

D. Numerical Results

The computation of an eco-driving cycle for HEV is done in the two steps described in Section V-B. The values of D f in Table IV are used for the two steps. For the first step, the DP mesh parameters are the same as in Section IV-B. For the energy management strategy computation, the following mesh parameters are used: ∆ξ = 0.025% for the SOC and a step of 2N.m for the engine torque.

The fuel consumption and the desired final SOC for scenario I are given in Table IX. Figure 12 summarizes the extra-fuel consumption in [%] with respect to scenario I. The various scenarios are compared with the same desired final SOC (ξ re f ) and with the same final time t f .

VI. CONCLUSION

The computation of eco-driving cycles for electric, conventional and hybrid electric vehicles using receding horizon control was addressed and formulated as an optimal control problem. The proposed method is based on dynamic programming. A first objective was to investigate the choice of the optimization frequency in order to find a trade-off between the optimality of the solution and the computation time of the algorithm. A relation between the optimization frequency and the look-ahead distance was identified and validated numerically.

Once D f had been optimized, a second objective was to study the impact of the look-ahead distance on the computation of eco-driving cycles. The result is that, for the three architectures under consideration, in the case of urban trips (mean speed values less than 40km/h), a look-ahead distance of 0.3 -0.5km leads to a sub-optimality in the energy saving less than 0.6% compared to the global solution. For high mean speed value trips, a look-ahead distance of 1 -1.5km can be used as the induced sub-optimality is less than 1%. Moreover, the computation times are relatively small: less than 1.2s for a conventional vehicle and less than 0.5s for the electric case. These results represent a first step toward a real-time implementation of the proposed algorithm in an Engine Control Unit (ECU).

Fig. 1 .

 1 Fig.1. Specific fuel consumption SFC (g/kWh) of the internal combustion engine as a function of the engine torque and the engine rotational speed. For confidentiality reasons, the data are normalized.

Fig. 2 .

 2 Fig. 2. Electric Machine Power [kW] as a function of the electric machine torque and the rotational speed. For confidentiality reasons, the data are normalized.

Fig. 3 .

 3 Fig. 3. Relation between β and the mean speed value v: Electric case.

Fig. 4 .

 4 Fig. 4. Predictive algorithm scheme for eco-driving.

Fig. 5 .

 5 Fig. 5. Sub-optimality and β for various values of D f [km] for D 0 = 2km.

Fig. 6 .

 6 Fig. 6. Extra-energy consumption [%] and the error on the final time [s] versus scenarios of Table IV: Electric case.

Fig. 8 .

 8 Fig. 8. Speed versus time for the EUDC: Conventional case.

Fig. 9 .

 9 Fig. 9. Corrected extra-fuel consumption [%] versus scenarios of Table IV: Conventional case.

Fig. 10 .

 10 Fig. 10. Speed versus distance for the A. Urban: Conventional case.

Fig. 11 .

 11 Fig. 11. Parallel mild-hybrid architecture.

Fig. 15 .

 15 Fig. 15. Predictive Eco-driving algorithm: inputs and outputs.

A sensitivity study using many driving cycles has been conducted to choose the best values of ∆v and ∆T e to find a trade-off between the solution optimality and its computation time: with the chosen parameters, the computation time is divided by 15 while the sub-optimality is lower than 0.5% comparing to the case with ∆v = 0.01m/s and ∆T e = 1N.m.

The extra-fuel consumption increases when the look ahead distance D 0 decreases, and its maximum value is less than 3.1% compared to scenario I (the sub-optimality remains in the same order of magnitude as the conventional case). Based on the nature of the trip, two cases are distinguished:

• Driving cycles with a low mean speed value (the ECE-15 and the Artemis Urban): the sub-optimality is less than 0.5% for all the scenarios considered. • Driving cycles with a high mean speed value (the Artemis Rural, the WLTC and the EUDC): scenario V for D 0 = 1km and D f = 0.26km, where the suboptimality is less than 2%, is a good choice. The vehicle speed and the SOC trajectories for the A. Rural driving cycle are given as a representative example in Figures 13 and14.

• The speed trajectories are very similar except for the case where D 0 = 0.5km. Small drops in the vehicle speed appear. These drops are related to the speed freeconstraint at the end of each segment. • The SOC trajectories are quite different. When the driving cycle is known in advance, the SOC variation is the largest (16% between the maximum and the minimum values of the SOC). When the look-ahead distance decreases, the maximum SOC variation decreases. This behavior is due to the constraint [START_REF] Van Keulen | Velocity trajectory optimization in hybrid electric trucks[END_REF] on the SOC at the end of each segment. To summarize, Figure 15 shows the inputs, the outputs and the parameters of the optimization algorithm (for all the architectures).

• β is estimated based on the mean speed value of the trip as explained in Section II-B.3 by equation [START_REF] Maamria | On the use of dynamic programming in eco-driving cycle computation for electric vehicles[END_REF]. The information about the mean speed can be provided by the GPS or any other source. The results are promising with an acceptable error on the final time of 4%.


Edwin Solano-Araque
email: edwin.e.solano-araque@renault.com

Guillaume Colin
email: guillaume.colin@univ-orleans.fr

Guy-Michel Cloarec
email: guy-michel.cloarec@renault.com

Ahmed Ketfi-Cherif
email: ahmed.ketfi-cherif@renault.com

Yann Chamaillard
email: yann.chamaillard@univ-orleans.fr

Determining vehicle acceleration from noisy non-uniformly sampled speed data for control purposes

Keywords: Vehicle acceleration, Controller Area Network, Non-uniformly sampled signal processing, Local Polynomial Regression, Savitzky-Golay Filter, Computational complexity

Vehicle acceleration is an important variable for many automotive control applications. In this paper, we present an approach to estimate vehicle acceleration from vehicle speed data coming from the CAN (Controller Area Network) bus. The proposed method, which can be seen as an extension of the Savitzky-Golay filter to non-uniformly sampled signals, re-samples them to a constant period while filtering noise coming from different sources and also provides a proper estimation of vehicle acceleration. We also consider the frequency response of the filtering effect of the method. Finally, some practical considerations for efficient implementation of the algorithm are given.

INTRODUCTION

Vehicle acceleration has a major influence on many automotive applications. For example, it is a determining factor in drivability [START_REF] Bellem | Objective metrics of comfort: developing a driving style for highly automated vehicles[END_REF]) and also in the energy consumption of electric vehicles [START_REF] Mruzek | Analysis of parameters influencing electric vehicle range[END_REF], [START_REF] Solano-Araque | Energy analysis of eco-driving maneuvers on electric vehicles[END_REF]). Therefore, acquiring a vehicle acceleration signal from real-driving field tests is of great interest for control design. These data can be obtained by installing acceleration sensors on the vehicle and obtaining this information directly. However, in some cases, this type of information may not be available; for example, when one wants to estimate the acceleration of the preceding vehicle for an automatic driving function, or when one must work in a given ECU of a car which has no access to the accelerometers but it has access to the vehicle speed. Another possibility is obtaining the vehicle speed signal, which is usually already present in the vehicle, and numerically deriving it in order to get an estimation of the acceleration. Nevertheless, given the fact that the vehicle's speed signal is usually very noisy, it can be hard to obtain an exploitable acceleration estimation from it.

In this work, we will consider the problem of getting an adequate acceleration estimation from the speed signal available on the CAN (Controller Area Network) bus of a commercial vehicle. This speed can correspond either to the controlled vehicle or the predecessor vehicle.

The CAN bus is widely used for automotive applications. This is due to its deterministic resolution of the contention (i.e. for a given information to be transmitted, the message to be sent is deterministically defined), low cost, and simple implementation [START_REF] Di Natale | Understanding and using the controller area network communication protocol: theory and practice[END_REF]). Today, most all (if not all) of the vehicles produced in Europe are equipped with at least one CAN bus [START_REF] Davis | Controller area network (can) schedulability analysis: Refuted, revisited and revised[END_REF]).

However, CAN characteristics lead to some phenomena that can pollute the signal being sent. One effect of this pollution is the fact that the sample time can vary slightly between two transmissions or, worse, some messages can be lost during transmission. This results in CAN signals not being uniformly sampled. This can cause problems when one wants to estimate the derivative of such a signal (e.g. the derivative of a speed signal).

In this paper, we propose a technique allowing us to get an exploitable estimation of vehicle acceleration, suitable to be used for control design. This technique is based on Local Polynomial Regression and can also be seen as an extension of the Savitzky-Golay filter [START_REF] Schafer | What is a savitzky-golay filter?[END_REF]) to non-uniformly sampled signals.

SIGNAL DISTORTIONS DUE TO CAN

PROTOCOL As already mentioned, the CAN protocol has some characteristics that may introduce noise in the signals. The main ones are :

-Priority based arbitration: The CAN arbitration protocol is priority-based, specifically it uses Carrier Sense Multiple Access/Collision Resolution (CSMA/CR) to determine access [START_REF] Davis | Controller area network (can) schedulability analysis: Refuted, revisited and revised[END_REF]). This means that if two nodes are trying to send a frame (sequence of bits) at the same time, the one with the highest priority (with a numerically lower identifier) will circulate and the other one has to wait for the next opportunity to be transmitted. Hence, it may introduce a non-deterministic delay in the signals being sent. -Bit stuffing: in order to avoid ambiguities in the frames being sent [START_REF] Davis | Controller area network (can) schedulability analysis: Refuted, revisited and revised[END_REF]), to allow nodes to synchronize their timing and to allow error detection [START_REF] Watterson | Controller area network (can) implementation guide[END_REF]), some bits can be added to the frames. This operation affects the length of the frames and, therefore, it may introduce a slight delay in the transmission. -Limited signal resolution: the maximum data length for a CAN frame is 8 bytes [START_REF] Di Natale | Understanding and using the controller area network communication protocol: theory and practice[END_REF]). This means that the maximum resolution of a signal being transmitted through a single frame is 64 bits. And in practice, the resolution is lower due to the fact that automakers usually send many signals packaged in one frame. -Loss of messages: On the CAN protocol there is no guarantee that the message that a node sends will be received by all of the others. Some messages may not be received by one or many of the concerned nodes. This would mean a loss of points in the transmitted signal. Hence, the method we propose allows a filtered uniformly sampled signal and a proper estimation of its derivative to be obtained from a noisy non-uniformly sampled signal.

METHOD FOR

The method consists of several stages, in particular: 1) defining the output signal time vector and "creating" timewindows, 2) getting the signal points associated to each window, 3) (local) polynomial approximation of window points, 4) reconstructing the point at t = t end [k] and estimating the signal derivative at this point.

Each of these phases is now presented in turn.

Output signal definition and windows consideration

First, it is necessary to define the time instants t end ∈ R R R Nsign 1 , for which we want to reconstruct the raw signal, whose time vector is t ini ∈ R R R Nsign 0 . In order to get a uniformly-sampled signal, the output-signal time points must be uniformly distributed, with a sampling time τ e . We will consider that the signal starts at t ini 1 = t 0 and it ends at t ini Nsign0 = t f .

Once the vector t end has been generated, we can consider the time-windows. First, we consider time-windows centered on each one of the points t end j and with a width τ win . Such an algorithm can be used to obtain information from an experimental setting a posteriori.

We will introduce the notation we are using. win j is the window corresponding to the i-th value of t end j . win j starts at t j win0 and ends at t j win f ; using equations 1 and 2 these values are calculated as a function of t end j and τ win , for a centered window.

t j win0 = t end j -τ win /2,
(1)

t j win f = t end j + τ win /2, . (2)
In order to promote that each window has enough points to apply the method, we also introduce the following constraints:

t end 1 t 0 + τ win 2 , (3)
t end Nsign1 t f - τ win 2 . ( 4
)
The following condition ensures a uniform sampling for the output signals:

t end j+1 -t end j = τ e , ∀j ∈ [1, ..., N sign1 -1]. (5)
By respecting these conditions, we can get an appropriate vector t end . For this paper, we have taken t end

1 = t 0 + τwin 2 and t end Nsign1 t f -τwin 2 . This is represented in Figure 1.

Obtaining the signal points associated to each window

From this point on, we are only working on the j-th window, win j . Each of the following stages can be applied to each of the windows.

At win j , we will get all the points between t win0 j t t win f j , each of which will be denoted x k j , k = 1, ..., m j . We will note x j the vector with all the points belonging to win j , as presented in (6).

x j = [x 1 j (t 1 ), ... , x k j (t k ), ... , x mj j (t mj )] T | ∃x k j (t k ) ∈ x j , ∀t k ∈ T j T j = {t| t ∈ t ini , t win0 j t t win f j } (6)

(Local) Polynomial approximation of the points within the window

Once we have gotten x j , it is possible to define a function f (t) to approximate the values of the points x k j (t k j ). For the method we propose we use a polynomial function, f np (t), as presented in (7), where n p is the order of the polynomial,

P np (t) = [1, t, t 2 , . . . , t np ] T and θ = [θ 0 , θ 1 , . . . θ np ] T . f (t) = np i=0 θ i .t i = P np (t) T θ (7)
In order to simplify the notation, from now on, we will omit the super-index j from most of the equations. f (θ, t) is a linear function w.r.t. θ, thus we can consider f (t)

w i n d o w s
Fig. 1. Time-windows for the proposed method as a linear measurement model [START_REF] Boyd | Convex optimization[END_REF]). Hence, by considering ν k as the error between the measurement x k and the estimation f (t k ), where t k ∈ t, we get:

x k (t k ) = P np (t k ) T . θ + ν k (8)
By considering all the m points in x, we find (9), where x ∈ R R R m are the points within the window,

T p = [P np (t 1 ), P np (t 2 ), . . . , P np (t m )] T ∈ R R R m×np is a ma- trix
where each row corresponds to the 0 th -n p th powers of t k , θ ∈ R R R np is a vector with the parameters of f (t), that we want to estimate, and ν ∈ R R R mj is the estimation error for each point. We will assume that m j n p .

x = T p . θ + ν (9)
By considering that a function of the form f (θ , t), such as the one given in 7, can "perfectly" (i.e. goodly enough for the application purposes) represent the signal behavior within the considered window, and that the difference between the measurements and the estimations corresponds to Gaussian noise (i.e. ν follows a Gaussian distribution), finding the parameters that best represent the points within the window becomes a least-squares problem [START_REF] Boyd | Convex optimization[END_REF]). Therefore, the optimal value for the parameters of f (t) for the window, θ , are given by (10).

θ = T T p . T p -1 .T T p .x (10)
3.4 Signal Reconstruction at t = t end j and signal derivative estimation To find the filtered signal value associated to win j , we need to find the value of f (θ , t) at t = t end j , this value will be associated to the output signal. This operation is presented in Figure 2.

Fig. 2. Output signal reconstruction

We can find an estimate of the signal derivative at t = t end j by applying (11):

dx dt (t end j ) = [0, 1, 2, . . . , n p ] diag (θ )         0 t end j 0 t end j 1 . . . t end j np-1         (11)
Once we have done so, we only have to apply the method to each of the windows to get a filtered signal and an estimate of its derivative.

Causal implementation of the method: "left-oriented" window

The method we have just presented, based on a centered window, is non-causal, which means that it is not possible to use it for a real-time implementation. However, one can easily get a causal version of this algorithm by using a window that only considers past values of the signal; it is equivalent to saying that t end j = t win f j . We may call such a window a left-oriented window. In that case, the values of t j win0 and t j win f can be calculated by applying equation ( 12) and the condition t

win f j = t end j . t j win0 = t end j -τ win . (12)
The conditions to promote having enough points to apply the method will become:

t end 1 t 0 + τ win , (13)
t end Nsign1 t f . ( 14
)
When one directly applies the polynomial approximation presented in section 3.3, the frequency response of the filter that one can get would not be acceptable for an application. This is due to the fact that a polynomial regression tends to be closer to the central values of the signal and yields greater errors for the corners of the approximation. However, a way to compensate for this fact is to weight the signal points by the absolute value of the difference between t end j , and each point time, t ini k . It is equivalent to minimizing the value of ν T W 1 ν, where

W 1 = diag (w), w 1 k = |t end j -t ini k |, which is a weighted least-squares problem.
In order to prioritize even more the points that are closer to the time at which one would like to reconstruct the signal, one can generalize the previous method by considering the magnitude of the l-th power of each gap. It is expressed by equation ( 15), which is a variant without loss of generality. The fact of dividing all the numbers by the l-th power of the maximum value of |t end j -t ini k | permits to normalize the values of w k to 1, in order to avoid numerical problems.

W l = diag (w) , w l k = |t end j -t ini k | l max r=1,...,mj |t end j -t ini r | l . ( 15
)
The solution of the weighted least-squares problem is given by:

θ = T T p . W l . T p -1 .T T p . W l .x (16)
4. FREQUENCY RESPONSE OF THE FILTERING AND CHOICE OF THE METHOD PARAMETERS An important characteristic of a filter is its frequency response. In this section we present the frequency response of the proposed methods, which can be viewed as interpolation filters, obtained numerically by applying it to sinusoidal signals of different frequencies and comparing the amplitude of the input and the output signals in order to get the gain. It is important to note that, given the fact that the proposed methods are not necessarily linear filters, it is possible that this representation does not give a comprehensive representation of the filter behavior.

In Figure 3, we compare the frequency response of both the proposed non-causal and causal methods with the response of a classic technique, namely a Butterworth filter. The parameters used for the method are: τ win = 200ms and n p = 1 for the non-causal method, and τ win = 850ms, n p = 1 and l = 2, which lead quite close to a cutoff frequency at 6dB of 2.997Hz (i.e. the frequency at which the output signal amplitude is attenuated by 50.12% w.r.t. the input signal). The Buttherworth filter has been designed such as to get the same cutoff frequency. We can observe that the non-causal method preserves more of the signal energy below the cutoff frequency and, at the same time, it attenuates the signal better for frequencies above the cutoff frequency. Concerning the causal method, we can see that for a certain range of frequencies there is a gain greater than 1, which means that those frequencies will be amplified by the filter. Nonetheless, the filter preserves (or amplifies) a great amount of signal energy in the range of interest, while better attenuating the noise higher frequencies.

PRACTICAL CONSIDERATIONS FOR OFFLINE

ALGORITHM IMPLEMENTATION When implementing the proposed method, the computational burden can be a problem. Probably, the main way of reducing it is by improving the way one solves the leastsquares problem for the polynomial approximation. However, another interesting way of improving the complexity is to reduce the burden associated to getting the points within each one of the time-windows. This is particularly true for long signals because, as is shown below, this is a task whose complexity scales quite rapidly as the input or output size increases. In this section, we present an appropriate algorithm to recover the points within the window while limiting the computational burden.

First we will consider a naive algorithm to get the point within the windows, presented in Algorithm 1. It will take as an input the original signal time vector, t, the number of inputs of the output signal, m, the sampling time of the output signal τ e and the size of the windows, τ win . The algorithm returns, for each window, a vector containing the indexes of the points within the window.

Algorithm 1 Determining the window points -naive version

Input:

Original signal : t ∈ R R R n ; Number of windows (output signal size) : m ; Method parameters : τe , τ win Output:

m vectors v j with the indexes of t belonging to the window j 1: begin Create m void vectors v j = ∅ , j = 1, ..., m 2:

for each value of the signal, i, do 3:

for time window, j, do 4:

if t[i] is contained by the window j then 5:

Append the value i at the end of v j 6: end if 7:

end for 8:

end for 9: end In order to analyze the complexity of the naive algorithm, we will use the pseudo-code of Algorithm 2, for which some programming choices have been done, including the use of dynamic memory allocation. The blue comments show the computational complexity of each task. To calculate the computational burden of a loop, one must add the cost of every task within the loop and, then, multiply it by the number of iterations in the loop, in "big O" notation.

For the naive algorithm, the main complexity is associated to the two nested loops, whose complexity is O(n.m), which means that the computational burden increases linearly with the input and the output signal sizes n and m, respectively. Given the fact that for a constant sampling time τ e , n and m are closely related, one can roughly consider it as a quadratic growth.

Concerning the memory cost of the algorithm, the algorithm requires generating m vectors, v j , with a mean size ṽ, roughly constant (in fact, ṽ depends on different parameters of the method, and it is indirectly associated

Algorithm 2 Determining the window points -naive version pseudo-code

Input:

Original signal : t ∈ R R R n ; Number of windows (output signal size) : m ; Method parameters : τe , τ win Output:

m vectors v j with the indexes of t belonging to the window j 1: begin 2:

Create m void vectors v j = ∅ , j = 1, ..., m O(m) 3: t win0 = 0; t winf = 0 O(1) 4:

for i := 1 : n do loop: O(n) 5:

for j := 1 : m do loop: O(m) 6:

t win0 := (j -1) * τe O(1) 7:

t winf := (j -1) * τe + τ win O(1) 8:

if t[i] t win0 and t[i] t winf then O(1) 9:

Append i at the end of v j O(1) (big) 10:

end if 11:

end for 12:

end for 13: end to m, but here we are neglecting this dependency). This represents a memory cost O(ṽ.m) = O(m).

The improved algorithm we propose is presented in Algorithm 3. It takes the same inputs as the naive algorithm, but it only returns two vectors of size m: v ini idx and v end idx , which respectively contain the initial and the final indexes of the points associated to each one of the windows.

Algorithm 3 Determining the window points -proposed algorithm

Input:

Original signal: t ∈ R R R n ; Number of windows (output signal size): m ; Method parameters: τe , τ win Output:

2 vectors with the initial and final indexes for each window, respectively v ini idx and v end idx 1: begin 2: Create and initialize the auxiliary scalars I ini , I end I ini 1 , t sign , t win0 , t win0 , t winf 3:

Create and initialize the vectors v ini idx and v end idx 4: for each signal value, i (starting at i = 2) do 5:

Save I ini 's last value in I ini 1 6:

Find the first window I ini whose t

I ini win 0 t[i]
7: Find the first window I end whose t

I end win f > t[i] 8: for k from Idx ini 1 + 1 to Idx ini do 9: if t[i] t win k
f and a point belonging to the window have not been found yet then 10:

i is the first index of the window k i.e. v ini idx [k] := i 11:

end if 12:

end for 13:

for k from I end to I ini do 14:

i is the last-until-now valid index of the signal within the window k i.e. v end idx [k] := i note: if in future iterations another index greater than i is found to belong to the window, v end idx [k] will be overwritten

15:

end for 16:

end for 17: end

In order to analyze the complexity of the proposed algorithm, we will use the pseudo-code of Algorithm 4, for which also some programming choices have been done; however, this time we have chosen to use static memory allocation, as we know the output size a priori.

Algorithm 4 Determining the window points -proposed algorithm pseudo-code

Input:

Original signal: t ∈ R R R n ; Number of windows (output signal size): m ; Method parameters: τe , τ win Output:

2 vectors with the initial and final indexes for each window, respectively v ini idx and v end idx 1: begin 2: As we can see on Algorithm 4, after the initialization phase there is a loop containing four independent internal loops. Concerning the computational burden, the initialization tasks (lines 2-4) can be executed in constant time O(1); the external loop is executed O(n) times, and each one of the internal loops is executed in a mean time O( m l ), l = 1, 2, where m l does not depend directly on the number of windows m. Therefore, by adding these results, we will get an amortized computational complexity

I ini = 1; I end = 1; I ini 1 = 0 O(1) 3: t sign = 0; t win0 = 0; t winf = 0 O(1) 4: v ini idx = [1 ; 0 m-1 ]; v end idx = 0 m O (
O(n(2 m 1 + 2 m 2 ) + 1) = O(n( m 1 + m 2 )), with ( m 1 + m 2 )
m when m is big.

Concerning the memory cost of the algorithm, we know exactly the amount of memory that the algorithm requires. We need six scalar variables (Idx ini , Idx end , Idx ini 1 , t sign , t win0 , t winf ) and two m-sized vectors (V idxini et V idxend ). Hence, the memory cost associated to the algorithm is O(2m + 6) = O(m). 4 presents the raw and processed signals of vehicle speed and acceleration. For the acceleration estimate, the following estimation methods were compared: applying Euler's method to the raw signal, applying Euler's method to the speed signal processed by a Buttherworth O1 filter, filtering with a Buttherworth filter the derivative of the raw signal, applying the proposed method with a centered window (non-causal method) and applying it with a left-oriented window (causal method). In order to fairly compare the methods, we have use the same parameters of the filters shown in Figure 3. In order to apply the linear filters, we have linearly interpolated the signals in order to get an uniformly-sampled signal (which introduces a distortion in the signal as it behaves as a low-pass filter).

We can notice that all the considered methods produce a much less noisy signal than the estimation based on the raw signal. However, the linear method introduces a delay in the speed signal.

Fig. 4. Vehicle Acceleration Estimation

Figure 5 presents a zoomed view of the velocity and accelerations. We can see that there is a real improvement in the acceleration due to the processing. However, all the causal methods (the proposed causal method and the linear ones) introduce a delay in the estimation; however, by using the left-oriented window, one gets a less noisy signal. Nonetheless, the proposed causal method also amplifies the signal values at some points. Also, there is a pronounced jump in the raw estimations at the beginning of the vehicle maneuver. This is due to the fact that commercial vehicle speed sensors are usually not able to produce a reliable estimation at low speed. Therefore, the speed signal has a value 0 when the vehicle speed is below a certain level. This produces a jump in the signal which has quite a strong effect on the acceleration estimation. We can notice that, when the non-causal proposed method is applied, this phenomenon practically disappears without introducing a delay in the signal.

CONCLUSIONS AND DISCUSSION

The proposed non-causal method provides an adequate estimate of vehicle acceleration without having an accelerometer, which could be used for many a posteriori Fig. 5. Vehicle Acceleration Estimation -Zoomed view applications. However, it is hard to say if the causal proposed method really improves the acceleration estimation; it could be the case when one needs to get a non-noisy rough estimate of the acceleration.

Also, an efficient way to implement a crucial step of the method has been proposed. This method could also be applied to some other problems where finding whether some given points belong to an equally spaced set of overlapping intervals is required.

In this work, we have considered a speed signal based on wheel sensors. Another possibility that could be applied to single-gear electric vehicles would be to estimate acceleration based on the electric machine speed; this is a more noisy signal but, at the same time, it does not have the problem of a low dead-zone and it could have a better resolution. The proposed method could properly filter such a signal in order to better exploit the information to get a relevant acceleration estimate.

Fig. 3 .

 3 Fig. 3. Frequency response of the proposed method and of a linear interpolation filter; (a) linear-scale, (b) logarithmic scale

Table . 1

 . presents a synthesis of the complexity results detailed above.

	Algorithm	Time	Memory
	Naive	O(n.m)	≈ O(m) (dynamic)
	Proposed	≈ O(n)	O(m) (static)

Table 1 .

 1 Comparison of naive and proposed algorithms in time and memory complexity 6. ESTIMATION OF VEHICLE ACCELERATION Finally, we applied a posteriori the proposed method to the vehicle speed signal obtained through sensors installed on each of the wheels, and transmitted by a CAN bus Figure

Orléans, France, June 23-27, 2019