Keywords: flow-structure interactions, instability, wakes

researchers and makes it freely available over the web where possible.

Three-dimensional mode selection of the flow past a rotating and inline oscillating cylinder David Lo Jacono 1,2, †, Rémi Bourguet 1 , Mark C. Thompson 2 and Justin S. Leontini 3

Introduction

When considering the flow-induced vibration of a circular cylinder in a free stream, vortex-induced vibration (VIV) is the only process capable of resulting in large-scale oscillations of the structure. During VIV, the body motion is predominantly in the cross-stream direction, even if the structure can respond in the streamwise direction, or rotate about its axis [START_REF] Jauvtis | The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[END_REF][START_REF] Horowitz | Vortex-induced vibration of a rising and falling cylinder[END_REF]. However, small departures from this basic flow t hat b reak t he r eflection symmetry, such as modifying the body to have an angle of attack, adding shear or turbulence to the incoming flow, o r f orcing t he c ylinder t o r otate ( as d one i n t his p aper) c an give rise to completely different bluff body flow-induced v ibration p henomena. Examples of such phenomena are transverse galloping and aeroelastic flutter [START_REF] Parkinson | The square prism as an aeroelastic oscillator[END_REF][START_REF] Zhao | Fluid-structure interaction of a square cylinder at different angles of attack[END_REF]. Such phenomena do not necessarily result in essentially crossstream oscillation, like VIV does, and the body motion can have a large streamwise component.

A particular example of such a streamwise flow-induced v ibration i s p resented in the recent study from [START_REF] Bourguet | In-line flow-induced vibrations of a rotating cylinder[END_REF] (BL15), which showed that a cylinder elastically mounted yet constrained to oscillate only in the streamwise direction could execute extremely large oscillations (with amplitudes over 2.5 times the cylinder diameter) if it was forced to constantly rotate. The change in symmetry of the boundary conditions, and the input of energy, excite flow-induced vibration modes that are not accessible in the non-rotating cylinder case. BL15 observed a number of exotic vortex configurations i n t he w ake o f t he b ody a nd a lso r ecorded a variety of three-dimensional structures.

Here, an attempt is made to explain the appearance of these various threedimensional structures using Floquet stability analysis. The underlying problem is simplified b y i mposing t he c ylinder r otation a nd a p urely s inusoidal streamwise oscillation, rather than allowing the oscillation to be driven by the flow. This simplification a ppears j ustified as th e co upled mo tion re corded by BL 15 is very close to sinusoidal. Various combinations of amplitude and frequency are run to encompass the range of parameters observed in the free case. Two three-dimensional modes are found, each with a distinct spanwise wavelength, and each linked to a unique instability mechanism. Depending on the frequency, either mode can be the first t o b ecome u nstable w ith i ncreasing a mplitude. C omparison w ith t he fully three-dimensional fluid-structure i nteraction s imulations o f B L15 s hows t hat t he first mode to become unstable sets the wavelength and subsequent saturated dynamics of the flow. T he f act t hat t he l eading m ode i s a f unction o f t he f requency e xplains why the three-dimensional wavelength varies as a complicated function of the rotation rate in the fully coupled simulations.

Section 2 describes the problem set-up and the numerical methods employed. Section 3 presents the results of the stability analysis and § 4 provides some concluding remarks.

Methodology

Base flows

The simulations were conducted using a well-validated spectral-element method [START_REF] Thompson | Three-dimensional instabilities in the wake of a circular cylinder[END_REF]. For the base flow, t he c ode s olved the incompressible Navier-Stokes equations in a frame of reference travelling with the centre of the cylinder

∂u ∂t = -(u • ∇)u - 1 ρ ∇p + µ∇ 2 u -ẍ, ∇ • u = 0,    (2.1)
where u is the velocity field, t i s t ime, ρ i s t he fl uid de nsity, p is th e pr essure field, µ is the fluid d ynamic v iscosity a nd x¨ = [ x ¨ 0 ] T i s t he a cceleration o f t he f rame of reference, equal to the linear acceleration of the cylinder linear motion given by

x = A d sin(2πf d t), (2.2)
where x is the displacement from the origin in the streamwise direction, A d is the amplitude and f d is the frequency of oscillation.

Since the circular cylinder is axisymmetric, the cylinder rotation was implemented by simply imposing the tangential velocity at the cylinder surface as a Dirichlet boundary condition. The boundary conditions for the velocity at the other boundaries were as follows. At the boundaries upstream and lateral of the cylinder, a Dirichlet boundary condition was used to impose the free stream velocity minus the frame velocity. At the downstream or outlet boundary, a Neumann boundary condition was imposed with the value of the gradient of the velocity in the normal direction set to zero.

For the pressure, a high-order Neumann condition was imposed at the body and the upstream and lateral boundaries, with the value of the pressure gradient in the normal direction derived from the Navier-Stokes equations [START_REF] Gresho | On pressure boundary conditions for the incompressible Navier-Stokes equations[END_REF][START_REF] Karniadakis | High-order splitting methods for the incompressible Navier-Stokes equations[END_REF]. At the downstream or outlet boundary, a Dirichlet condition was imposed, with the value of the pressure set to zero.

The spectral-element method solved equations (2.1) in the weak or variational form, similar to all finite-element-based methods. The domain was decomposed into 916 quadrilateral elements, that were free to have either straight or constant-radius curved sides. These elements were further internally decomposed by using seventh-order Lagrange polynomials as shape functions. These Lagrange polynomials were associated with Gauss-Legendre-Lobatto quadrature points for efficient calculation of the spatial integrals required for the weak form solution. A thorough description of spectral-element methods can be found in [START_REF] Karniadakis | Spectral/hp Element Methods for Computation Fluid Dynamics[END_REF].

Temporal discretization employed a three-way time splitting scheme (e.g. see [START_REF] Thompson | Hydrodynamics of a particle impact on a wall[END_REF]. Here, the integration from the start to the end of a time step was split into three substeps: the first substep integrated the advection and frame acceleration terms (the first and fourth terms on the right-hand side of (2.1)) to a velocity field at an intermediate time using an explicit second-order Adams-Bashforth scheme; the second substep integrated the pressure term from this first intermediate velocity field to a second intermediate time using an implicit second-order Crank-Nicolson scheme; the third substep integrated the diffusion term from this second intermediate field to the velocity field at the end of the step using an implicit second-order Crank-Nicolson scheme.

To solve for the pressure, the divergence of the second substep equation was taken, and the incompressibility constraint was imposed on the velocity field at the end of the substep. This results in a Poisson equation that can be solved for the pressure field, which can then be used to complete the second substep.

Details of this spectral-element scheme can be found in the textbook from [START_REF] Karniadakis | Spectral/hp Element Methods for Computation Fluid Dynamics[END_REF].

Previous studies have shown that the domain size can be important in flows such as those studied here. The domain size has been fixed throughout the study at Ω ∈ [-20D : 90D][-55D : 55D], similar to other recent studies of similar flows [START_REF] Rao | Three-dimensionality in the wake of a rapidly rotating cylinder in uniform flow[END_REF][START_REF] Bourguet | In-line flow-induced vibrations of a rotating cylinder[END_REF]. The mesh resolution has been assessed by varying the order of the polynomial shape functions of the spectral elements. The results show that base flow quantities vary by less than 0.1 %, and eigenvalues from the Floquet stability analysis vary by less than 1 % when varying the polynomial order from seven to eight.

Perturbation fields a nd s tability analysis

The Floquet methodology relies on a periodic base flow. A ll t he simulations undertaken here happen to be periodic within the parameters tested. The first s tep for the global Floquet stability analysis is to form equations that govern the evolution of small perturbations. This is done by decomposing the velocity and pressure fields into base flow a nd p erturbation c omponents, s ubstituting t his i nto t he e quations o f motion outlined in (2.1), subtracting the terms for the base flow, a nd l inearizing (discounting nonlinear terms). The resulting equations are

∂u ∂t = -((u • ∇)u + (u • ∇)u) - 1 ρ ∇p + µ∇ 2 u , ∇ • u = 0,    (2.3)
where u is the perturbation velocity field and p is the perturbation pressure field. This system can be further decomposed by considering perturbations that are harmonic in the spanwise direction. Doing so results in a set of equations that decouple for different spanwise wavelengths, λ, and the wavelength appears explicitly in the equations [START_REF] Barkley | Three-dimensional Floquet stability analysis of the wake of a circular cylinder[END_REF]. These equations were solved using the same spectral-element method as outlined above. Dirichlet boundary conditions stating that the perturbation was zero were imposed at the cylinder surface and the upstream and lateral boundaries. The zero-normal-gradient condition was imposed at the outflow boundary. For the pressure perturbation, the condition for the normal pressure gradient was the same as that imposed on the base flow. The analysis progressed by first evolving the base flow to a periodic state, then the base flow and perturbation flow were integrated in time together.

Floquet stability analysis progresses by looking for eigenvectors and associated eigenvalues of a linear operator L that maps a system from one period to the next. In the context of the flow considered here

u (t + T) = Lu (t), (2.4)
where T is the period. In implementation, the operator L is never explicitly formed; its action is calculated by integrating the equations governing the perturbation evolution (2.3) forward by one period. Further details of performing linear stability using this time-stepper approach can be found in [START_REF] Tuckerman | Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems[END_REF] and details of the implementation used here are presented in [START_REF] Leontini | Stability analysis of the elliptic cylinder wake[END_REF]. Since the operator is not explicitly formed, its corresponding eigenvectors and eigenvalues cannot be found directly. However, the leading eigenvectors (those with the largest eigenvalues) can be found through iteratively applying the operator, i.e. integrating over multiple periods. Here, the leading eigenvectors and eigenvalues are found using Arnoldi iteration.

These leading eigenvalues are the most relevant from a stability point of view. Equation (2.4) shows that the eigenvalues can be interpreted as the ratio of the perturbation field f rom o ne p eriod t o t he n ext. I f t here i s a n e igenvalue ( here referred to as the Floquet multiplier, µ f ) such that |µ f | > 1, then the perturbation field is growing from one period to the next, and the base flow i s u nstable. H ence, t he aim of the analysis is to determine for which parameters the base flow i s s table (all |µ f | 1) or unstable (at least one |µ f | > 1). The eigenvector associated with a given eigenvalue is often referred to as a Floquet mode -it is basically the perturbation flow field.

Problem set-up

The basic problem investigated is that of a circular cylinder immersed in a free stream. The cylinder is rotating at a constant rate about its axis, while simultaneously performing harmonic oscillations in the streamwise direction. Simulations were conducted for flows w ith c onstant R eynolds n umber Re = UD/ν = 1 00 a nd normalized rotation rate α = ωD/U = 3, where U is the free stream velocity, D is the cylinder diameter, ν is the kinematic viscosity and ω is the rotation rate. The normalized amplitude A d = A/D and frequency f d = fD/U, where A and f are the amplitude and frequency of streamwise oscillation respectively, were systematically varied over the ranges 0 < A d < 2.4 and 0.04 f d 0.11.

Results

Floquet multiplier spectra

Figure 1 presents the magnitude of the Floquet multiplier, |µ f |, as a function of the spanwise wavelength, λ. Each panel shows data for a single frequency of oscillation, while each line represents data for a constant amplitude. It is clear that for the lower frequencies, there are two clear peaks in these Floquet multiplier spectra. Each peak can be associated with a distinct mode of instability. One of these modes (that focused around the shorter wavelengths and is marked with red lines and points in the figure) has a real, positive Floquet multiplier. The other (marked with blue lines and points in the figure) h as a r eal, n egative F loquet m ultiplier, i ndicating i t i s a subharmonic mode which repeats over two oscillations of the body.

As the frequency is increased, the peaks in the spectra for both modes shift towards shorter wavelengths. Typically (except for the case at f d = 0.06) the subharmonic mode becomes unstable first, i .e. h as a F loquet m ultiplier |µ f | > 1 a t a l ower a mplitude than the shorter-wavelength mode. However, |µ f | for the shorter-wavelength mode increases very rapidly with amplitude, and if the amplitude is increased beyond the critical point where |µ f | = 1, the value of |µ f | for the shorter-wavelength mode quickly exceeds that of the subharmonic mode.

Spatial structure and mechanism of instability of the two modes

Example instantaneous images of each of the modes are presented in figure 2 . The contour lines in the figures r epresent vorticity o f t he b ase fl ow an d th e co lour contours represent levels of the spanwise vorticity of the perturbation field.

For the shorter-wavelength harmonic mode, the perturbation is focused in the layer of fluid d irectly a djacent t o t he b ody. T ogether w ith t he c losed s treamlines a round the body, this is taken as evidence that a centrifugal instability mechanism plays a role in the destabilization of this mode.

For the longer-wavelength subharmonic mode, a potential instability mechanism is not as clear. There appears to be some growth very close to the body, again suggesting a centrifugal mechanism. However, the subharmonic nature of the mode suggests some interaction between wake vortices of the same sign separated in time by one period [START_REF] Yildirim | Mode C flow transition behind a circular cylinder with a near-wake wire disturbance[END_REF], generated by the break in symmetry of the wake stemming from the body rotation. The low-frequency example image shown in figure 2 s upports t his i dea, w ith g rowth o f t he p erturbation i n the wake vortices. However, this evidence of growth in the wake vortices is not as clear at higher frequencies. While it seems likely the same cooperative mechanism (i.e. centrifugal instability near the body supplemented by further growth via interaction of wake vortices) exists at the higher frequencies, a definitive c onclusion c annot be drawn. Lines and points marked in red represent data for the harmonic mode, those in blue represent data for the subharmonic mode. Filled black points mark a common amplitude for which one of the modes is close to marginally stable. Red points mark example cases shown in figure 2 . ( a-c) S how t hat t he s ubharmonic mode is first q uenched, t hen a mplified wi th in creasing fr equency. Ha rmonic mo de da ta is not presented for the higher frequencies in (d-f ) as the mode is not resolved in the range of wavelengths tested at these frequencies.

Critical values of amplitude and wavelength

Figure 3 plots the critical amplitude (the lowest amplitude at which |µ f | = 1) and the critical wavelength (the wavelength for which |µ f | = 1) for both the harmonic and subharmonic modes as a function of the oscillation frequency f d . Figure 3(a) shows that typically, the subharmonic mode becomes unstable at a lower amplitude than the (a) show that typically the subharmonic mode becomes unstable at a lower amplitude of oscillation; except in the range 0.05 < f d < 0.07, where the harmonic mode becomes unstable first with increasing amplitude. The data in (b) show that the wavelength of the harmonic mode is consistently shorter than the subharmonic mode. harmonic mode, implying that for amplitudes of oscillation at least up to the critical amplitude, |µ f | for the subharmonic mode is greater than |µ f | for the harmonic mode. The exception to this is in the range 0.05 < f d < 0.07, where the harmonic mode is the first to become unstable with increasing amplitude. For comparison with the freely oscillating case, this is important as it might be expected that as the amplitude of the free case develops from zero, the first mode to become unstable will be the mode that will go on to the fully saturated three-dimensional flow. Of course, this is not universal -nonlinear effects that will occur as the amplitude moves further from the critical amplitude could mean the fully saturated flow bears little or no resemblance to the linear mode calculated here. However, if there is a bifurcation sequence leading to the fully saturated flow, the first step in this sequence will be governed by which of the linear modes is unstable first w ith i ncreasing a mplitude, s ubject t o e qual projection of the background noise onto the modes.

Figure 3 (b) shows that the critical wavelengths of the two modes are typically well separated, and both decrease with increasing frequency. The separation of the wavelengths of the two modes presents an interesting prediction: if there is a situation in the freely oscillating case where the frequency can vary, it may be possible to switch the mode which leads the transition, and this may result in a change in the wavelength of the flow. T he c omparison b elow i n § 3 .4 s hows t hat t his i s i ndeed the case.

Around f d = 0.06, the subharmonic mode is quenched. This is due to the fact that as the amplitude is increased beyond the critical value for the harmonic mode, the growth rate or Floquet multiplier for the harmonic mode increases very rapidly -much more rapidly than the Floquet multiplier for the subharmonic mode. This high sensitivity to the amplitude also supports the idea that this mode is driven by the centrifugal instability of the layer of fluid a djacent t o t he c ylinder. I magine a perturbation introduced as a wave packet at a location inside the layer of fluid inside the closed streamlines that encircle the body. For the perturbation to grow, this wave packet must grow over the time it takes to make one orbit around the cylinder [START_REF] Bayly | Three-dimensional centrifugal-type instabilities in inviscid two-dimensional flows[END_REF]. Therefore, the time scale for the multiplication of this growth rate is the orbit time, which will be similar to the rotation time of the cylinder itself. For the cases tested here, α = 3. A frequency of f d = 0.06 gives a period of T = 16.66, implying the cylinder makes 50 rotations per oscillation cycle. A small increase in the ratio of the perturbation from one rotation to the next therefore results in a huge increase in the ratio of the perturbation from one oscillation to the next.

Comparison with the freely vibrating case

Here, these controlled oscillation stability results are used to interpret threedimensional DNS simulation results of a rotating cylinder that is elastically mounted and free to oscillate in the streamwise direction presented by BL15. In this previous study, it was reported that the clear three-dimensional structures were formed, and that there was no clear relationship between the spanwise wavelength of these structures and the frequency, amplitude, or reduced velocity U * -U/( f n D), where

f n = √ k/m/(2π)
, k is the spring stiffness and m is the cylinder mass. For four values of U * = 9, 11, 15.5 and 23, the amplitude and wavelength measured during the freely oscillating simulations are plotted as a function of the measured response frequency in figure 3 . T he r otation p arameter α = 3 i s t he s ame i n t he freely oscillating cases and the stability cases presented here. It is clear that the amplitude of the response is just enough to exceed the critical amplitude required for threedimensional instability. It is also clear that the measured wavelength is reasonably well predicted by the critical wavelength of the first m ode t o b ecome u nstable. A striking feature is the large drop in wavelength measured in the DNS between the cases at f d = 0.04 and f d = 0.06 (U * = 23 and U * = 15.5). This corresponds to a change in the first m ode t o b ecome u nstable f rom t he s ubharmonic m ode t o t he h armonic mode.

This change in mode is also highlighted in figure 4 . S pace-time d iagrams a re shown for the DNS simulations, plotting contours of the spanwise velocity as a function of time and spanwise distance z, along a line located at (x, y) = (5, 0) for the same four values of U * as presented in figure 3 . T hese d iagrams s how t hat i n t he D NS, t he only case that repeats on the oscillation period (1-periodic) is that at U * = 15.5 -all the others repeat over two cycles, i.e. all the others are subharmonic (2-periodic). Therefore, the lack of any regular relationship between the amplitude, frequency and reduced velocity of the freely oscillating case modelled in the DNS simulations is explained. The wavelength selected in the freely oscillating case is a function of which mode is the first to become unstable, and this is a function of frequency of oscillation.

3.5. Relationship to the non-oscillating rotating cylinder [START_REF] Rao | Three-dimensionality in the wake of a rapidly rotating cylinder in uniform flow[END_REF] presented an in-depth study of the three-dimensional modes present in the wake of a rotating cylinder that is not oscillating over a range of rotation rates covering the value α = 3 used in the present study. Due to the low frequencies used here, a comparison with the limiting case of no oscillation is warranted.

However, the two situations -oscillating and non-oscillating -have major differences in flow topology; the primary one being that the oscillating flows here exhibit periodic vortex shedding, and the non-oscillating cases do not. In fact, the non-oscillating cases do not exhibit a traditional wake at all, as the high-speed rotation effectively wraps all the vorticity produced at the cylinder surface around the cylinder, and the flow r emains s teady. T hese d ifferences r ule o ut a q uasisteady explanation of the oscillating case. Regardless, two modes grow on the steady flow: m ode E , a steady-steady bifurcation, with a critical spanwise wavelength of approximately 1.3 for Re = 100; mode F, a centrifugal instability with a critical wavelength less than 1.

The similarity in wavelength between mode E and the subharmonic mode found here may suggest a similar instability mechanism, which [START_REF] Rao | Three-dimensionality in the wake of a rapidly rotating cylinder in uniform flow[END_REF] conjectures may be a combination of centrifugal and hyperbolic effects. This is plausible in the oscillating case too, and § 3.2 discusses the potential for a cooperative mechanism.

The centrifugal mode F instability is also very similar to the harmonic mode found here. As this mode is an instability of the region of flow e nclosed b y t he streamlines orbiting the cylinder, it is almost unaffected by the vortex shedding in the wake, and hence its appearance in both the oscillating and non-oscillating case is not surprising. [START_REF] Rao | Three-dimensionality in the wake of a rapidly rotating cylinder in uniform flow[END_REF] also presents further evidence in this regard, showing that the trend of critical Re for mode F is practically unaffected by the onset of a low-frequency vortex shedding mode at higher rotation rates. Therefore, while not quantitatively similar, it appears that the instability mechanisms that occur in the non-oscillating case are qualitatively similar to those found here for the oscillating case.

Concluding remarks

Floquet stability analysis of a rotating cylinder performing streamwise oscillations in a free stream has been conducted. Two distinct modes, with different wavelengths and distinct mechanisms of instability, have been identified, and both of them can lead the transition to three-dimensional flow. The stability results are used to interpret the measurements of three-dimensional structures in three-dimensional DNS simulations of a rotating cylinder that is freely oscillating in the streamwise direction -the stability results are able to explain that the non-monotonic relationship between the response amplitude, frequency and reduced velocity is due to a change in the linear mode which leads the transition to three-dimensional flow, which is a function of the frequency of oscillation.

FIGURE 1 .

 1 FIGURE 1. Magnitude of the Floquet multiplier |µ f | as a function of spanwise wavelength λ. Each curve represents a unique value of the amplitude of oscillation A d . Each panel contains data for one frequency of oscillation f d : (a) f d = 0.05; (b)f d = 0.06; (c) f d = 0.07; (d) f d = 0.08; (e) f d = 0.09; ( f ) f d = 0.11. Lines and points marked in red represent data for the harmonic mode, those in blue represent data for the subharmonic mode. Filled black points mark a common amplitude for which one of the modes is close to marginally stable. Red points mark example cases shown in figure2. ( a-c) S how t hat t he s ubharmonic mode is first q uenched, t hen a mplified wi th in creasing fr equency. Ha rmonic mo de da ta is not presented for the higher frequencies in (d-f ) as the mode is not resolved in the range of wavelengths tested at these frequencies.

FIGURE 2 .

 2 FIGURE 2. Spatial structure of the two modes. Red/blue colour contours represent positive/negative spanwise vorticity of the perturbation field while solid/dashed contour lines represent positive/negative base flow vorticity. (a) Short-wavelength harmonic mode for wavelength λ = 1.0. (b) Longer-wavelength subharmonic mode for wavelength λ = 1.6. Amplitude A d = 1.5 and frequency f d = 0.05 in both cases.

  FIGURE 4. Space-time diagrams taken by recording the spanwise velocity at (x, y) = (5, 0) over time from the three-dimensional DNS simulations for the four sample cases. Red/blue contours represent positive/negative velocity. The solid line at z = 4 is simply to highlight the behaviour over time at a given spatial location. Only the case at U * = 15.5 repeats every period, showing the change from the subharmonic mode to the harmonic mode. (b) Two pairs of images of isosurfaces of ω x at one instant, first for U * = 11 and then for U * = 15.5. Each pair consists of an image taken from the three-dimensional DNS simulations and an image reconstructed from the stability calculations at similar amplitude and frequency. Arrows indicate correspondence between the space-time diagram and the isosurface image. These images show the structure of the flow is well predicted by the linear Floquet mode.

† Email address for correspondence: david.lojacono@imft.fr c

Acknowledgements

This work was financially supported by the Australian Research Council Discovery Projects scheme via grant DP150103177. Part of this work was performed using HPC resources from CALMIP (grants 2017-P1248 and 2018-P1248).


Philippe Fabian
email: philippe.fabian@univ-eiffel.fr

Abderrezak Rachedi

| Cédric Guéguen

Selection of relays based on the classification of mobility-type and localized network metrics in the Internet of Vehicles

Keywords: Internet of Vehicles, Machine Learning, Selection of Relay, Vehicle-to-Vehicle

The ever increasing amount of connected mobile devices and users in the Internet of Vehicles (IoV) and the fact bicycles, electric scooters and users' smartphones are also connected poses a real challenge in terms of ensuring Quality of Service (QoS).

The constantly changing topology due to the high-mobility of devices and users greatly impacts its stability and the connectivity of devices. Furthermore, to ensure safety of people in the case of autonomous vehicles it is of paramount importance to ensure excellent reliability. This is why we have developed a solution that is based on Machine Learning (ML) to classify devices according to their mobility profile and that uses a scoring system to select the best candidates to act as Mobile Relays (MR) amongst devices with a suitable mobility profile. The scoring system allows to find critical locations in terms of user density. This solution does not require a dedicated infrastructure such as Road Side Units (RSUs). Simulations results will show the proposed solution increases the packet delivery ratio by up to 6%, reduces the energy consumption by up to 30% and increases the efficiency of bandwidth usage without sacrificing the end delay of users and devices compared to the state of the art.

INTRODUCTION

Mobile connected devices are more and more present 1 2 3 : pedestrians with mobile phones, motorized vehicles equipped with an on board computer or public connected bicycles and scooters. All these devices form the Internet of Vehicles (IoV). They can potentially connect to one another or they can directly connect to Internet using for instance either a 4/5G connection to a cellular antenna or a WiFi connection through an Access Point (AP). However, in such a dynamic network of mobile devices it is challenging to ensure topology stability and thus Quality of Service (QoS). Indeed, the topology is constantly changing because users and devices travel in different directions at different velocities, they connect and disconnect from the network according to their needs or to save battery life, or because objects such as vehicles cause multipath fading, resulting in some devices gaining or losing connectivity. Furthermore, users and devices have different needs including high throughput for a multimedia stream on a smartphone, minimum delay for autonomous vehicles to avoid a collision or low energy consumption to ensure enough remaining residual energy on a device until the next recharge.

It is not always possible for a device to have a direct connection to a cellular antenna because one might not be available from a given Internet Service Provider (ISP) in low populated areas (deployment costs too high), the user has no subscription plan (too expensive or because the user is traveling abroad) or because the radio conditions are bad (path loss, shadowing, multipath fading). It is possible to resolve all those issues by using relays to forward data from and to an end user/device. Several solutions propose the deployment of static relays (dedicated architecture such as Road Side Units (RSU)), Mobile Relays (MR) 4 5 or the usage of parked vehicles [START_REF] Rachedi | BadZak: An hybrid architecture based on virtual backbone and software defined network for Internet of vehicles[END_REF] . Static relays have limitations such as deployment cost, lack of flexibility and mandatory regular maintenance of the equipment to ensure reliability. On the other hand, MRs allow for more flexibility and can lower deployment costs. In our vision, a dynamic solution using regular devices as relays is the next step compared to dedicated MRs as it allows to implicitly adapt the number of relays and select new ones whenever radio conditions or the needs of users and devices evolve and has the advantage of no maintenance costs. This solution also offers great flexibility by removing the need for a dedicated infrastructure. It allows devices without cellular access to connect to Internet through others. Local communications need not to be forwarded to cellular antennae, which can potentially reduce overhead. Furthermore, this solution can be deployed on top of an existing routing protocol 7 8 9 for even better QoS.

This is why we propose in this paper a solution that dynamically selects relays amongst devices present in the area of interest based on a Machine Learning (ML) classification algorithm. The proposed solution is an extension of our previous work [START_REF] Fabian | Dynamic selection of relays based on classification of mobility profile in a highly mobile context[END_REF] where the selection of a dynamic relay is not only based on the classification of mobility but also on localized metrics such as the proposed Expected Packet Count (EPX) and average delay of the different locations. This makes the new solution even better because taking into account localized network metrics allows to select relays in the most critical locations in order to better service users and devices. Furthermore, compared to our previous work, we expanded the types of users in our simulations by adding pedestrians. The highlights of this contribution are :

• We extended the classification algorithm compared to our previous work [START_REF] Fabian | Dynamic selection of relays based on classification of mobility profile in a highly mobile context[END_REF] by adding localized network metrics computed after the classification algorithm has been run in order to take into account the most critical locations in real-time when selecting relays;

• We created a new localized metric, the EPX which represents the packet density, in order to find critical areas where the density of devices is highest;

• The performance evaluation has been further improved by adding pedestrians (along with bicycles and cars) in scenarios considering a light and high density of users, allowing us to test the solution in a very realistic way.

The remaining of this paper is organized as follows : the state of the art is presented in section 2, the proposed solution is explained in section 3, section 4 details the performance evaluation and, finally, the paper is concluded in section 5.

RELATED WORK

Existing solutions presented in this section are regrouped into three categories : general, ML based and Cluster Head (CH) based. General approaches (subsection 2.1) concern the selection of relays based on a mathematical model. They are based on Markov chains or cost functions. ML based approaches (subsection 2.2) use machine learning to select relays. Finally, CH based solutions work by electing one relay for each group (or cluster) of users, hence the appellation "cluster head" (subsection 2.3).

General approaches

Authors in [START_REF] Abada | An Adaptive Vehicular Relay and Gateway Selection Scheme for Connecting VANETs to Internet via 4G LTE Cellular Network[END_REF] propose a solution to select mobile relays amongst vehicles to allow other vehicles to connect to Internet through them. Neighboring vehicles are connected through WiFi and relays connect to Internet using LTE. Several metrics are used to select relays such as speed, Received Signal Strength (RSS), stability of links and distance between vehicles. Though this solution is good, it does not consider several mobility types and simulations are done in straight lines on highways. In [START_REF] Liu | Opportunistic mobility utilization in flying ad-hoc networks: A dynamic matching approach[END_REF] , the authors use Matching Game to predict the radio conditions of Unmanned Aerial Vehicles (UAVs). Once UAVs have determined their own needs, the proposed solution dynamically adapts the mode selection, the time scheduling and the channel allocation to improve performance on the network. However, the authors of this paper studied cases where the number of UAVs is no more than a few tens and only a subset of drones generate data packets. This solution is thus not adapted to a road network where several hundreds of users can generate data at the same time. The solution proposed in [START_REF] Dao | Resource-aware relay selection for inter-cell interference avoidance in 5G heterogeneous network for Internet of Things systems[END_REF] uses a relay selection mechanism to help nodes that suffer from high interference. If the radio conditions of a node are too bad, it will scan its surroundings to find suitable neighbors to relay its data. The best neighbor is found using a M/M/1/K first-come-first-serve scheme. This solution allows any node to potentially serve as a relay but it is studied in a context without mobility, which does not make it adapted to the IoV. The authors in [START_REF] Liao | Mobile relay deployment in multihop relay networks[END_REF] propose to add MRs to help fixed relays during heavy traffic periods. Markov chains are used to predict when heavy traffic periods will happen. With this information, the optimal paths for MRs are computed and they are sent on those paths to assist fixed relays. Although this solution considers mobility of users, MRs are considered to move at a constant speed and end users can only connect to a relay (fixed or mobile). This limits the flexibility of the proposed solution as any user must be in range of a relay and there is a required dedicated architecture (fixed relays).

Machine Learning based solutions

Mchergui et al. [START_REF] Mchergui | Relay Selection Based on Deep Learning for Broadcasting in VANET[END_REF] propose the use of deep learning to select relays in order to increase performance. The solution is studied on a road network in the context of Vehicular Ad hoc Network (VANET). RSUs are used to allow vehicles to connect to Internet. The combination of several features in a deep learning algorithm is interesting. However, RSUs are static and this makes the proposed solution non-adapted in an environment where it is not always possible to deploy a dedicated infrastructure. In [START_REF] Zhou | A fuzzy-rule based data delivery scheme in VANETs with intelligent speed prediction and relay selection[END_REF] authors use a ML model based upon K-nearest neighbors for relay selection in the context of VANETs. The purpose is for a vehicle to select the best RSU in range or, if none are nearby, to select the best candidate (i.e. another vehicle) to forward data to a relay. Though the use of several features and multi-hop data forwarding makes the model adapted to VANETs, a dedicated architecture of RSUs is required. Furthermore, only one class of mobility (motorized vehicles) is considered. Authors in [START_REF] Dhurandher | Reinforcement Learning-Based Routing Protocol for Opportunistic Networks[END_REF] propose a Fuzzybased Q-Learning Routing Protocol to find the best next hop based on the node's residual energy, movement and buffer space. The proposed solution is tested in different scenarios with several tens of nodes with velocities up to 1.5 m/s. However, in the context of IoV there can be up to hundreds of users and they can move significantly faster. Authors in [START_REF] Abdelreheem | Deep Learning-Based Relay Selection In D2D Millimeter Wave Communications[END_REF] use deep learning for the selection of a relay in the case of a node having bad radio conditions due to shadowing. There are several metrics used to determine the best relay. However, this solution only considers a user can connect to a base station through one relay (one hop).

Electing a relay through cluster head selection

The solution proposed in [START_REF] Zhang | Data Forwarding at Intersections in Urban Bus Ad Hoc Networks[END_REF] exploits the proximity of buses at intersections to forward data. Bus density, road connectivity and the path of bus lines are used to determine on which bus the data is forwarded. This solution relies on the presence of a bus network with a minimum density to efficiently forward data. Also, the delivery delay from simulations ranges from 7 to 10 seconds. These make the proposed solution inadequate in areas without bus lines or for applications needing a real-time latency. Al-Kharasani et al. [START_REF] Zukarnain | An Adaptive Relay Selection Scheme for Enhancing Network Stability in VANETs[END_REF] propose a clustering solution to connect nodes in a VANET. Several metrics like bandwidth, connectivity, velocity and distance are used to create a cluster of vehicles. However, the proposed solution is tested in scenarios using constant speed and only motorized vehicles (i.e. high-mobility users) are considered. In [START_REF] Zhao | A scheme of D2D-based delay analysis and vehicle relay algorithm in VANET[END_REF] , the authors propose a scheme where vehicles with similar motion are regrouped in clusters. The cluster head is also connected to the infrastructure using a cellular link. The clustering allows to regroup vehicles with a similar trajectory together so only the head need to be connected to the infrastructure or to other clusters' head. However, in the context of urban mobility several different types of users are present and many intersections will cause users to travel apart from one another with a high probability. In [START_REF] Baidas | Joint relay assignment and energy-efficiency maximization in energy-harvesting downlink/uplink clustered nonorthogonal multiple-access networks[END_REF] , authors tackle the problem of joint relay assignment and energy-efficiency maximization to assign a relay to clusters of users. To do so, authors propose to formulate the problem considering users' constraints as a maximization problem, solve the problem for each pair of cluster/relay, and finally assign the relay. The assignment is achieved with a linear complexity. However, the performance evaluation focuses on relatively static users and only 12 of them are considered. This is non-adapted to the IoV where hundreds of devices with very different mobility patterns are present at the same time. In [START_REF] Dahmane | WeiSTARS: A weighted trust-aware relay selection scheme for VANET[END_REF] , authors propose a scheme to select relays based on network and trust metrics in VANETs to select trustworthy relays. However, the proposed solution makes use of dedicated RSUs and only one class of mobility is considered.

All of these solutions either lack support for several mobility types typically found in cities (pedestrians, e-scooters, bicyles and motorized vehicles), require a dedicated architecture to operate or they consider a limited amount of hops to connect devices to base stations. In our vision, a more general approach with a greater flexibility is mandatory in improving system performance. Thus, we propose a solution that supports different mobility types, can adapt to devices' needs in real-time and only requires one cellular antenna to work. Few devices with cellular connectivity are needed to serve as relays.

DYNAMIC SELECTION OF RELAYS

In this section, we present the architecture and an overview of the proposed solution (sec. 3.1), we discuss the MR selection algorithm in details (sec. 3.2), we analyze selected algorithms (sec. 3.3) and we explain the novel metric EPX (sec. 3.4).

Architecture and overview

The proposed architecture, inspired from Software Defined Networking (SDN) [START_REF] Taleb | On multi-access edge computing: A survey of the emerging 5G network edge cloud architecture and orchestration[END_REF] , is depicted on figure 1 and it is made up of two parts : the Fog part and the Cloud part. Users with Internet connectivity (pedestrians with a smartphone, connected cars and so on) and the cellular antenna (all circled in green on fig. 1) compose the Fog part of the network. On figure 1, fog users have a wireless link to the antenna represented as a black-dashed line. Internet and the Control Server are the Cloud part of the architecture (in light blue on fig. 1). This architecture allows a higher degree of flexibility because a dedicated infrastructure is not required, except for one cellular antenna. Even without an antenna users can still communicate locally. Given we assume all of them can be mobile at any time, such an architecture is well adapted compared to a static one where devices can potentially leave the range of the antenna.

The algorithm will run either using a two step mode (initialization mode) or a four step mode (on-going mode). The initialization mode is as follows :

1. Users that have Internet access connect to the controller to send their mobility data (such as velocity, acceleration, type of road the device is on, etc.) as depicted on figure 1; 2. The controller uses the data to determine the mobility profile of the users. The ones with the appropriate profile are selected as relays.

This mode is executed when the network is yet to be built. This is similar to our previous work [START_REF] Fabian | Dynamic selection of relays based on classification of mobility profile in a highly mobile context[END_REF] . Once the network is built, the selection mechanism can be further optimized by running the on-going mode :

1. Users send to the controller their mobility data (such as velocity, acceleration, type of road the user is on, etc.) and their network metrics (such as the average delay of their packets and the average amount of transmitted packets), along with their location as depicted on figure 1; 2. The controller uses the mobility data to determine the mobility profile of users. Users with the appropriate profile are pre-selected and tagged as potential relays;

3. The controller uses the network metrics to find the most critical locations. A critical location is one where the values of metrics are considered as too high; 4. Relays are elected on priority in critical locations. Then, as long as there are locations considered as critical other relays are selected as well. Figure 2 is an example of the topology once the algorithm has been run and relays have been elected. Relays are users directly connected to the antenna. We remind the reader that in this case there were four users with Internet access (see figure 1), labeled from 1 to 4 on figure 2. Here, bicycle user 1 and pedestrian user 2 are the elected MRs. Pedestrian user 3 was not elected because her radio conditions are too bad (she is located inside a building). Car user 4 was not elected as a MR because he is a highmobility user, not suitable for the stability of the topology. The "appropriate profile" is determined according to local needs. In our case, we consider soft-mobility users (i.e. pedestrians) have the best profile because of their low mobility. Indeed, if they are elected as relays the resulting topology is expected to be more stable compared to electing faster users like connected cars.

The addition of localized metrics such as average delay and EPX (see subsection 3.4) allows to find locations where the density of users and devices is higher. This allows to select relays at these locations which minimizes the number of hops needed by users and devices to reach the relay. The probability of re-transmissions is lower and the end-to-end delay is likely to decrease. This in turn increases QoS. Limiting the number of elected relays minimizes interference between them. Please note that the localized metrics are not part of the classification. They are used after the running of the ML algorithm to enhance the selection mechanism by providing the critical locations.

In our case, the features to classify the mobility-type of users are the position, the velocity and the maximum acceleration. As explained earlier in this section, this data is sent to the controller. However, only the position through time is required (velocity and acceleration can be deduced from it). We suppose the data is measured on a device each 100 ms and sent every second to the controller. This can be further optimized by sending only one position value as long as it does not change. One x, y, z position can fit in three 4-octets variables (so 10 ⋅ 12 for 1 second). This data needs to be sent once every tens of seconds in practice because the selection process happens at about this frequency, so the overhead is not very high (a few packets per minute).

Mobility classification and optimal relay selection

The proposed solution is made up of two principal parts :

1. The ML algorithm classifies users according to their mobility profile and 2. Their network metrics allow to determine which ones from the above step are the most appropriate to act as relays.

We define a network composed of users. Each user possess one vector of mobility-related metrics :

= [ 1 , 2 , … , ], (1)
where is the number of mobility-related metrics. These metrics can be the current or maximum velocity of the user, his location (i.e. the type of road he is on) and so on. The user periodically records these metrics. Each user also stores one vector of network-related metrics :

= [ 1 , 2 , … , ], (2)
where is the number of network-related metrics the user is experiencing (such as delay, throughput, ETX and so on). The metrics are periodically recorded as well. Mobility-related metrics (eq. 1) are used by the ML algorithm as features for the classification. The accuracy in determining the mobility type of the user (e.g. pedestrian, bicycle, car, etc.) depends on the choice of features and specific algorithm used. We thus introduce the vector of features for each user :

= [ 1 , 2 , … , ] s.t. | | > 0 ∀ and ≤ . ( 3
)
That is, the features in equation 3 are a subset of the available mobility-related metrics from equation 1. The same features are used from each user and at least one feature must be used. These features are recorded for a length of time equals to . The purpose of is to choose the amount of time during which mobility-related metrics are gathered. Indeed, gathering metrics for a longer amount of time tends to increase accuracy, but it means waiting more before establishing the profile of the user. This can lower QoS, especially in time-sensitive applications so a compromise must be found between accuracy and the delay before predicting the class of a user. The mobility class predicted by the algorithm for each user is thus :

= ( , ). (4)
In equation 4, is the selected machine learning algorithm, are the features of user (see eq. 3) and ∈ is the predicted class of user . Given the prediction is not always correct, we introduce the real class of user , ∈ , such that :

= if the prediction is correct. ∈ ( ⧵ ) if the prediction is wrong. ( 5
)
If the prediction is correct, the predicted and real class of user are the same. On the other hand, a failed prediction will classify the user in the wrong class which is any other class than the real one. The best mobility class is determined as the most appropriate class to serve as relay. For instance, a slow-mobility class (e.g. pedestrian) can be desirable to increase topology stability.

Once the mobility class of users is determined, the network-related metrics (from eq. 2) are used to find the best location for electing a relay. First, we divide the studied area into per tiles. The dimensions of each tile can be as small as to have one user maximum per tile (e.g. 1m x 1m) or as big as to have only one tile for the whole area and for all users :

, = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1,1 1,2 ⋯ 1, 2,1 2,2 ⋯ 2, ⋮ ⋮ ⋱ ⋮ ,1 ,2 ⋯ , ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 6
)
where each , is one tile. We now define a scoring function to compute the score of each tile given the values of networkrelated metrics of users :

, = ( ) ∀ ∈ , . (7)
The purpose of tiles is to identify the most critical locations (e.g. tiles with high delay, high ETX, low throughput, etc.) to elect one or more relays inside them. The vector of the most critical locations is introduced :

⊆ , s.t. | | > 0. ( 8
)
That is, the vector of the most critical locations (eq. 8) is a subset of the matrix of scores (eq. 6). A location can be considered critical when QoS requirements are not met, for instance. Finally, relays can be elected once these locations are determined :

∈ iff ( ∈ , ) ∧ ( ∈ ) ∧ ( ∈ ) | , ∈ . ( 9
)
That is, a user will serve as relay if it is located in any critical cell (see eq. 8), has Internet access (belongs to ) and is of the correct (best) mobility class (eq. 4 and 5). Note that Internet access is essential for a node to be a potential relay because it must be able to forward traffic back and forth from/to Internet. That is, a node without Internet access is not even a sub-optimal relay compared to one with Internet access; it can not be a relay at all. Hence, these nodes are not considered by the ML algorithm. Note that in the case where a critical tile has no user that can serve as relay, these users will eventually be connected to a neighboring tile. This neighboring tile will become critical and relays will be elected to address the issue.

Algorithm 1

The selection relay algorithm. end if 19: end for 20: return Algorithm 1 shows the relay selection process. It is run whenever the maximum number of relays is not reached, for instance when a MR is disconnecting due to low battery life or because it does not need to be connected to the network any longer. The parameter of the algorithm (arg[0]) is the maximum number of relays. Lines 1 to 3 are the initialization steps. The gathering of users' metrics is done in lines 4 to 6 (eq. 1 and 2). In lines 7 to 11, the classification algorithm is run for each users to determine if its mobility profile is suitable to serve as a MR (eq. 4 and 5). If so, it is added to the set of potential relays (line 9). Then, in lines 12 through 19, the score of the potential relays is computed to find the ones that will best increase performance of the network (eq. 8). The score is computed in line 13 (eq. 7) and, if the user is satisfying, it is added to the set of relays in line 14 (eq. 9). The adding of relays can be interrupted if enough relays are present (line 16). Finally, the set of relays is returned in line 20. Note that a user must be in the correct location to become a potential relay. An important advantage of using ML is the fact it can select a relay amongst non-optimal candidates (for instance a car) if there are no optimal candidates (pedestrian) available. This can yield to good performance if the car has a temporary low mobility, perhaps due to a traffic jam. This allows more flexibility compared to an exhaustive search solution as we look at the current features (velocity, acceleration) of nodes rather than their overall mobility profile.

Benchmark of different ML classification algorithms

We tested several ML algorithms to classify the mobility-type of devices. The results are shown in table 1. The training set was generated with the Simulator of Urban MObility (SUMO, section 4.1). The size of the testing set is 40% of the total data [START_REF] Sakiz | A survey of attacks and detection mechanisms on intelligent transportation systems: VANETs and IoV[END_REF] .

The accuracy of the algorithms depends on the delay allowed to gather data. When the delay is very short (2 s for instance), the error is rather high because the user was not active long enough. For example, a car just starting might stop at a red light right away and it is not possible to differentiate it from a bicycle. However, although a longer delay allows for better accuracy, it also means a longer time before selecting a new relay. During this time, no connectivity will degrade QoS and might not be the right choice. This is why a compromise must be found between accuracy and delay to yield the best QoS possible.

Table 1 shows various values of accuracy for three different algorithms: K-nearest neighbors, Support Vector Machine (SVM) and decision tree. For each algorithm, two lines show the accuracy using a 10-cross-validation in classifying users as pedestrian, bicycle or car for the training and testing sets. Decision tree yields the best accuracy with 99% for the training set and ranging from 90% to 100% during testing. K-nearest is the second most accurate with a training accuracy of 97% and a testing accuracy from 80% (6 s) to 90% (14 s). SVM is the least accurate with correct classification in 93% of cases for the training set and between 73% and 85% when testing. The longest delays do not significantly improve the accuracy of the classification. This means a very long delay is not necessarily better. This is why we use decision tree with delays up to 6 s in the evaluation of performance (sec. 4). In terms of complexity, decision tree is restricted to a depth of 2 to avoid over-fitting and there are 3 labels as output (see section 4).

Table 2 depicts the precision, recall and the resulting f-score for the decision tree algorithm according to the type of mobility and delay allowed for the classification. The number of pedestrians, bicycles and cars ("Number of samples") used during the evaluation is also in the table. The proportion of bicycles is lower compared to the two others (less significant). Cars are easier to identify because of the potential high velocity they can reach, although when the allowed delay is very short (2 s), they do not have time to reach such speed. This explains the relative low recall (65%) compared to a classification delay of 14 s (recall is then 97%). The precision when identifying pedestrians is also quite low when the delay is 2 s (22%). This is coherent with results from section 4, where using a delay of 4 s or 6 s yields better results (the classification is more relevant).

New metric : the Expected Packet Count (EPX)

The EPX is different from Expected Transmission Count (ETX) [START_REF] Couto | High-throughput routing for multi-hop wireless networks[END_REF] . It is the total amount of (re-)transmissions needed for packets sent from a location to reach their destination during a certain time frame. Each hop in the path of these packets counts towards EPX including re-transmissions. It is the "density of packets" of a given location. A high EPX is the consequence of long paths to reach a relay, bad radio conditions (many re-transmissions) and/or simply high traffic.

A B C

P 1 EPX P1 = 3 (c) retries to send 1 to . EPX is incremented to 3.

A B C

P 1 EPX P1 = 3 (d) The packet reaches its final destination, . The packet reads an EPX of 3 and originates from .

FIGURE 3

An example of the computation of EPX: packet 1 is emitted from and has destination sink . The EPX, stored in 1 , concerns the area where is located.

Figure 3 shows a simple example of the computation of EPX. One packet is emitted from node and its destination is the sink node . First, when the packet is about to be sent after its creation, the value of EPX is incremented to 1 (fig. 3a). The location of the packet's source is kept in the packet. When the packet is received on , its EPX value is incremented as it is about to be forwarded to (fig. 3b). The transmission fails due to interference and is not acknowledged by . increments the EPX value again and resends the packet (fig. 3c) to . The value of EPX is now 3. Finally, successfully receives the packet (fig. 3d). Note that the value of EPX from 1 concerns the area where the source (here, ) is located, not the source itself. Thus, if a sink receives multiple packets from a same location, the total EPX is the sum of all EPX values from all packets from this area.

PERFORMANCE EVALUATION

The evaluation of performance is done by comparing four different algorithms. Two scenarios are studied and two parameters vary (thus, there are 4 series of graphs). The two first graphs of each series (figures a and b in section 4.3) represent the scenario with a low density of users (column Low density scenario in table 3). The third and fourth graphs of each series (figures c and d in section 4.3) refer to the scenario with a high density of users (column High density scenario in table 3). The line charts depict the variation of traffic generated by users, whereas the bar charts represent the change of selection delay (the time during which network related metrics are gathered for the purpose of classification of the mobility type of users). The studied algorithms are : Given the classification requires some time to gather data to determine the mobility profile of users, the selection delay of the algorithms without classification is always 0 s (no delay). The score function (see eq. 9) for the solution using classification with metrics combines the delay of tiles with the EPX.

Simulation setup

Simulations were done using the simulator Omnet++ 5.6.1. Users' trip were generated using SUMO (Simulation of Urban MObility) 1.3.1. The studied area is the location of the University Gustave-Eiffel (France) and was downloaded from OpenStreetMap. SUMO allows to randomly generate trips for different kinds of users (we used cars, bicycles and pedestrians) in a realistic way. Users follow roads from the downloaded OpenStreeMap area (sidewalks in the case of pedestrians) and obey to traffic rules. Once the trips of users have been created in SUMO, they are exported to a XML file. Then, the Omnet simulation will create one "module" (or node) for each user that was generated by SUMO. To train the ML algorithm doing the classification, we ran another simulation on a different location, allowing us to test the generalization of the algorithm. The delay of the prediction allows to gather more or less data from users. A lower delay is less accurate (about 70% accuracy for 2 s) because fewer data from a user's position, acceleration, speed and so on were accumulated. On the other hand, waiting a longer time to make the prediction means a higher precision (about 94% for 6 s) but users are not connected during that period. The list of the most important parameters is presented in table 3. The size of the area is about 1.96 km 2 . Users are created on the network each 3 s (low density case) or 0.75 s (high density case). The total and instantaneous number of users represent, respectively, the amount of users created on one run and how much of them there are at any time. About 5% of the instantaneous number of users are elected as relays. Most of the users are in cars and the smallest part of users are bicycles. Simulations were 1000 s in length with 10 repetitions to increase results accuracy. All users except relays generate data following a Variable Bit Rate (VBR). Data packets are set to the minimum size of an IPv6 packet [START_REF] Deering | Internet Protocol, Version 6 (IPv6) Specification[END_REF] and users have a maximum throughput of about 10 Mbps [START_REF] Stallings | Data and Computer Communications[END_REF] . The maximum range of users is about 105-110 m with a path loss based on a -10 ( ) ( is the distance) model 28 29 . The amount of static RSUs (fourth algorithm) is set to 5. This is about the same amount compared to the others (5% of relays), though it is purposely constant to reflect the fact the infrastructure is static.

Studied performance indices

This section provides explanations concerning the Performance Indices (PIs) used for the performance evaluation. We use the term "PI" to differentiate results shown in this section from metrics used to build the topology (as explained in section 3). PIs are presented in the same order as the next section (sec. 4.3) for the sake of clarity.

Packet delivery ratio

The packet delivery ratio (PDR) represents the percentage of generated data packets that reach their final destination. Some factors such as frequent disconnections of nodes from their neighbors or very bad radio conditions causing the amount of retransmissions exceeding the maximum number of retries will cause PDR to decrease. Note that a packet being re-transmitted many times on each hop can still yield to a high PDR, though other PIs such as delay will increase. The higher the PDR the better.

End delay

The end delay is the time taken for a packet to reach its final destination. It depends on how good radio conditions are (many retransmissions will increase delay) and on the amount of congestion (buffer occupancy). A low throughput means nodes require more time to send packets, increasing delay. The amount of nodes a packet has to go through on its path (hop count) does not necessarily increase delay, as a long path composed of nodes with high throughput and low congestion will perform better than a short congested path. The lower the delay the better.

Energy consumption

The energy consumption per node represents the amount of energy spent by each node on transmitting and receiving packets, including signaling and data packets. A high number of transmissions per unit time implies increased duty-cycling. Bad radio conditions resulting in many re-transmissions and topology instability causing more signaling are examples of causes of increase energy consumption. Energy consumption can be very important, especially for devices that can not have their battery recharged. The lower the consumption the better.

Throughput usage

The throughput usage represents the percentage of the total throughput a node can achieve that is in use. For instance, if a node can achieve a maximum throughput of 10 Mbps, a throughput usage of 20% means the node sends on average 2 Mbps of packets, including signaling and data. A lower value is achieved when good radio conditions (fewer re-transmissions) are present or if the topology is more stable (less signaling). The lower the throughput usage the better.

Relative amount of packets reaching destination

PDR, end delay and throughput usage are not enough to see how well the different solutions perform. Indeed, if a solution generates a very low amount of packets due to lack of connectivity, the PDR of these packets can still be very high, their delay very short and the usage of throughput will be low as well. The relative amount of packets reaching destination is the number of packets that reached the final destination compared to the total amount of packets created. It is the PDR times the total amount of generated packets and represents how well nodes are connected to a relay. The higher the better.

Simulation results

The following results are regrouped by the studied metric. We remind the reader that in line charts (where traffic load varies), the selection delay is 2 s. In bar charts (representing the variation of selection delay), the traffic load is set to 15 packets/s.

Packet delivery ratio

Figure 4 shows the PDR of the different algorithm. When studying the impact of traffic in the low density scenario, (fig. 4a), the proposed solution performs significantly better at higher loads. This is because taking into account network metrics allows it to adapt the location of relays where the need is greater. On the other hand, at lower traffic loads (figs. 4a) all solutions perform similarly, except the static RSUs. This is because mobile users eventually travel out of range of RSUs resulting in a very low PDR. In the case of a higher density of users (figs. 4c and4d), the proposed solution almost always performs significantly better because of more traffic generated in absolute. Indeed, in the case of a higher density of traffic, the use of network metrics allows to select relays in the most critical locations. Algorithm with classification has similar results compare to without classification on figures 4a and 4c because the low amount of relays (5%, see table 3) and the fact that users within the correct class of mobility are not differentiated from one another makes it harder to select a candidate in a critical location. The same is seen on figures 4b and 4d. The variance of the proposed solution increases with the selection delay in the low density scenario (fig. 4b) because the network metrics used for the scoring function are older (less accurate). This is not the case in the second scenario 4d because the higher density of users greatly increases the amount of network metrics and better locations are available to select a relay.

FIGURE 5

The average end delay results of the different scenarios. The vertical axis also has the same scale on all graphs.

End delay

Figure 5 shows the average end delay of data packets for the different solutions. It is interesting to see that although the PDR was not better in the low density scenario (fig. 4a) for the proposed solution, the end delay is better (fig. 5a). This is because even though there are fewer choices to select a potential relay, the network metrics still help in determining the most critical locations. When studying the selection delay in the low density scenario (see fig. 5b), the proposed solution is better, except when the selection delay is 4 s. As stated at the beginning of section 4.1, this is due to the fact an average selection delay will decrease the accuracy of the network metrics but it will only increase the classification accuracy by a moderate amount. On the other hand, either a short selection delay (accurate network metrics and shorter waiting time) or a longer selection delay (much more accurate classification) are preferable. For the high density scenario (figs. 5c and 5d), the delay is comparable to the other algorithms. In this case, the high amount of users makes the delay and EPX increase in all locations. The other solutions perform similarly (figs. 5a, 5b, 5c and 5d) for the same reasons mentioned before. The variance of static RSUs tends to be very high because the location of RSUs varies from one run to another and users tend to travel out of their range. This highlights the relevance of dynamically selecting relays in the context of IoV.

Energy consumption

The average energy consumption per user is depicted on figure 6. In all scenarios (figs. 6a, 6b, 6c and 6d), the proposed solution consumes less energy. The higher PDR (fig. 4) and the similar or lower end delay (fig. 5) show that there are less re-transmissions, resulting in a lower consumption of energy. It is worth noting that even though there are more users and packets generated on the high density scenario (see fig. 8), the consumption per user is not necessarily higher. This is due to the fact that more paths are

FIGURE 7

The average throughput usage in the different scenarios.

available, so the PDR is higher (as depicted on fig. 4c) and there are less re-transmissions. Static RSUs consume more energy because of lower PDR (fig. 4) and higher delays (fig. 5) caused by more re-transmissions and losses.

Throughput usage

Figure 7 shows the average throughput usage per node of the different solutions. In the low density and high density scenarios (figs. 7a and 7c), the usage increases with load for all algorithms. The usage is higher in the low density case compared to the high density case due to a higher PDR (fig. 4a compared to fig. 4c) and more re-transmissions. Static RSUs' PDR is similar in low and high density cases, so the throughput usage is similar. The selection delay (figs. 7b and 7d) does not impact much the throughput usage, though the variance increases slightly for the algorithm with classification and metrics in the low density case (fig. 7b) because of the decreased accuracy of the network metrics.

Relative amount of packets reaching destination

The relative amount of packets reaching their final destination compared to the amount of generated packets is presented in figure 8. The proposed solution shows better results in the low density scenario at higher loads (fig. 8a) and in the high density scenario (fig. 8c). This follows the results on the PDR (fig. 4a and fig. 4c). The selection delay does not impact the relative amount of arrived packets for the proposed solution in both scenarios (figs. 8c and 8d). The low PDR of Static RSUs explain why the amount of packets reaching the destination is low.

FIGURE 8

The relative amount of packets that reached their final destination.

CONCLUSION

Guaranteeing QoS in the IoV is still a challenging issue. The different mobility profile of users (pedestrians, bicycles, e-scooters, cars, buses and so on), their different velocity and their different paths contribute to the instability of the topology. We are convinced our solution to dynamically select relays amongst users by running a machine learning algorithm to classify users according to their mobility profile and the use of network metrics to find critical locations is key in improving QoS in the IoV.

Results show the proposed solution increases the PDR by up to 6%, reduces the average energy consumption of users by up to 30% and increase topology stability by reducing re-transmissions. Future work will focus on extending the solution to support more wireless technologies and network metrics.

FIGURE 1

 1 FIGURE 1 The proposed architecture made up of the Fog part and the Cloud part. Here, 4 users have access to Internet through the cellular antenna.

FIGURE 2

 2 FIGURE 2The established topology after a run of the proposed solution.

  As is about to transmit 1 , the ETX is incremented. received the packet and tries to send it to , thus incrementing the EPX. The transmission fails.

1 . 2 . 3 . 4 .

 1234 Without classification does not use any classification. Relays are selected randomly. In line charts it is represented by red circles. In bar charts it is represented as red-crossed bars with a selection delay of 0 s; With classification uses classification of mobility profile to select a new relay 10 . Relays with low mobility are favored. In graphs, this algorithm is represented by blue diamonds (line) and blue horizontally-lined bars (bar charts); With classification and metrics is the proposed solution. It uses mobility classification (as With classification) and a score function (see below). It is represented in line charts as green squares and in bar charts as green-filled bars; Static RSUs does not use classification (selection delay of 0 s). RSUs are placed on the studied as static units. In graphs, this algorithm is represented by purple crosses (line charts) and purple vertically-lined bars (bar charts).

FIGURE 4

 4 FIGURE 4The packet delivery ratio of the different scenarios. The vertical axis has the same scale for all 4 graphs.

( a )

 a The impact of the average network load in the low density scenario. The impact of the selection delay in the low density scenario. The impact of the selection delay in the high density scenario.

  The impact of the selection delay in the low density scenario. The impact of the selection delay in the high density scenario.

TABLE 1

 1 The precision in classifying the correct type of mobility against the allowed delay. The table shows the values for the training and testing sets using 10-cross-validation.

	Classification delay (in s)
	2	4	6	10 14

TABLE 2

 2 The precision, recall, f-score and the sample size when using the decision tree algorithm.

			Classification delay (in s)
			2	4	6	10	14
		Pedestrian 22	42	70	72	73
	Precision (in %)	Bicycle	100 100 100 100 100
		Car	100 100 100 100 100
		Pedestrian 100 100 100 100 100
	Recall (in %)	Bicycle	100 40 100 78 100
		Car	65	88	93	95	97
		Pedestrian 36	59	82	84	84
	F-score (in %)	Bicycle	100 57 100 88 100
		Car	79	94	96	98	96
		Pedestrian	2	5	7	13	19
	Number of samples	Bicycle	1	5	5	9	9
		Car	20	33	40	66	95

TABLE 3

 3 Values for the most relevant parameters. † The variable parameter is either traffic density (from 10 to 15 packets/s) or the relay selection delay (from 2 to 6 s). The solutions without classification always use 0 s).

	Parameter Value	Low density scenario High density scenario
	Playground size	1400 m ⋅ 1400 m	1400 m ⋅ 1400 m
	User generation period	≈ 3 s	≈ 0.75 s
	Num. of users, total	≈ 330	≈ 1250
	Num. of users, inst.	30-50	100-150
	Ratio of relays	max. 5 %	max. 5 %
	Car:bike:pedestrian ratio	3 : 1 : 1.25	2.75 : 1 : 1.5
	Simulation duration	1000 s	1000 s
	Repetitions	10	10
	Traffic model	VBR	VBR
	Traffic density †	10 to 15 packets/s	10 to 15 packets/s
	Relay selection delay †	[0, 2, 4, 6] s	[0, 2, 4, 6] s
	Data packet size	1280 bytes	1280 bytes
	Max. user throughput	10 Mbps	10 Mbps

  The average energy consumption per node in the different scenarios.

		3500							3500						3500							3000
				Without classification											Without classification
	Energy consumption per node (mW)	500 1000 1500 2000 2500 3000		With classification With classi. and metrics Static RSUs		Energy consumption per node (mW)	500 1000 1500 2000 2500 3000			Without classification With classification With classi. and metrics	Energy consumption per node (mW)	500 1000 1500 2000 2500 3000		With classification With classi. and metrics Static RSUs		Energy consumption per node (mW)	500 1000 1500 2000 2500			Without classification With classification With classi. and metrics
												Static RSUs													Static RSUs
		0	10	11	12	13	14	15	0	0	2	4	6		0	10	11	12	13	14	15	0	0	2	4	6
			Average traffic generation per node (packets/s)			Sink selection delay (s)				Average traffic generation per node (packets/s)			Sink selection delay (s)
		(a) The impact of the average	(b) The impact of the selection		(c) The impact of the average	(d) The impact of the selection
		network load in the low density	delay in the low density scenario.		network load in the high density	delay in the high density
		scenario.												scenario.						scenario.
	Average throughput usage per node (%)	22,5 25 27,5 30 32,5 35				With classification With classi. and metrics FIGURE 6 Without classification
						Static RSUs
		20	10	11	12	13	14	15
			Average traffic generation per node (packets/s)

We trained the algorithms with different ratios for testing/training sizes and the difference is marginal.