Julien Pansiot
email: julien.pansiot@inria.fr

Edmond Boyer
email: edmond.boyer@inria.fr

3D Imaging from Video and Planar Radiography

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The ability to capture intrinsic body structure in motion is of interest in a number of fields related to medical imaging such as computer-assisted surgery, biomechanics, and sports science. Many applications consider video or depth cameras and infer skeletal motion from surface observations using prior models. However, this strategy does not provide real measures on the internal structure and the estimated skeleton does not match the actual bone structure due to multiple factors such as inaccurate skeletal model and complex elastic tissue motion. With the aim to provide better measures to observe intrinsic structures in motion, and validation purposes, we investigate in this paper a new strategy that recovers dense 3D volumetric models of moving samples.

To this purpose, we combine a video-based surface motion capture system that provides motion cues, with a single static planar X-ray imaging device that captures the inner structure. We present the first step towards three-dimensional volumetric motion capture by investigating first rigidly moving samples, assuming limited prior knowledge on the captured samples. A key concept of our approach compared to traditional tomography is that it does not consider motion as low-amplitude noise to be corrected, but at the contrary as a source of information, ensuring the capture of X-ray images from multiple viewpoints.

As a result, the proposed configuration can consider moving samples as well as several sensors (eg. two X-ray devices). Yet less accurate than a CT-scanner, it yields a less expensive low-dose solution, taking benefit of equipment widely available in clinical environments. Our volumetric reconstruction method builds on super-resolution techniques [START_REF] Liu | On Bayesian adaptive video super resolution[END_REF] to optimally exploit X-ray samples and infer 3D attenuation. It relies on an X-ray image formation model (see Fig. 1) accounting for 2D sensor noise as well as 3D geometric errors. This model is associated with a volumetric L 1 smoothness prior to constrain the reconstruction, allowing for a limited number of input views. All these elements are integrated within a Bayesian framework for backward inference.

To summarize, the key contribution introduced in this paper is a Bayesian approach to 3D imaging of a moving sample which accounts for both sensor and calibration inaccuracies using a generative model.

Related Work

Currently, the two well-established classes of methods to recover attenuation models are planar radiography and Computed Tomography (CT). The former can capture motion at competitive frame rates with a low dose, but is limited to integrated two-dimensional information. The latter, in the form of Multi-Detector CT (MDCT) and Electron Beam CT (EBCT), can capture accurate dense 3D volumetric attenuation models at up to 30 fps, yet exposing patients to higher ionising radiations, present much higher costs, and limited versatility.

For these reasons, few-views cone-beam computed tomography (CBCT) techniques [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] have gained interest. Several methods have been devised to reconstruct relatively accurate 3D attenuation models from a limited number of cone-beam images [START_REF] Bang | CT reconstruction from a limited number of X-ray projections[END_REF][START_REF] Yang | A nonlinear least squares method for solving the joint reconstruction and registration problem in digital breast tomosynthesis[END_REF]. Nevertheless these methods are limited to static samples.

The combination of motion capture and emission tomography has also been investigated but either requires markers [START_REF] Mcnamara | A flexible multicamera visual-tracking system for detecting and correcting motion-induced artifacts in cardiac SPECT slices[END_REF], or is designed for low-amplitude motion correction [START_REF] Hutton | A hybrid 3-D reconstruction/registration algorithm for correction of head motion in emission tomography[END_REF]. Similarly, motion correction in tomography has been largely covered in traditional CT/PET [START_REF] Dawood | Lung motion correction on respiratory gated 3-D PET/CT images[END_REF] as well as CBCT [START_REF] Zhang | Correction of motion artifacts in cone-beam CT using a patientspecific respiratory motion model[END_REF]. Again, our strategy differs since we consider motion as a mean to vary viewpoints which helps capturing moving samples.

A method to reconstruct 3D attenuation from a limited number of arbitrary X-ray views was proposed in [START_REF] Sidky | Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[END_REF], but assumes reasonably good calibration. The general approach by [START_REF] Pansiot | Combined visible and X-ray 3D imaging[END_REF] reconstructs volumetric attenuation from a limited number of X-ray views of a moving object which motion is estimated using videos. While building on a similar concept, we take however a different and more formal generative strategy inspired by image super resolution [START_REF] Liu | On Bayesian adaptive video super resolution[END_REF]. In this domain, Bayesian inference methods have demonstrated their ability to optimally exploit noisy observations with uncertainty modelling and we extend them to our X-ray imaging configuration.

Generative image model

As mentioned, our generative model builds on existing image formation model [START_REF] Liu | On Bayesian adaptive video super resolution[END_REF] to explain the X-ray images given the 3D model. Our method takes as input a set of X-ray images of a rigidly moving sample. The images are first registered in a common framework using the motion estimated by a multi-view motion capture system. A dense attenuation model of the moving sample, represented as a voxel grid, is then reconstructed using the entire X-ray image sequence.

We detail below the main components of this model. In order to account for the multiple sources of noise present in the acquisition process, we introduce a generative image formation model, as illustrated in Figure 1. This model is associated with a sparse prior, ie. a TV-norm.

Image formation

We discretise the continuous absorbance problem in 3D as a weighted sum over the voxels v j along the given ray ω, d j being the distance covered within the voxel v j and µ j the attenuation assumed uniform within v j , defining the absorbance I(ω) in function of the emitted and transmitted intensities L 0 and L(ω):

I(ω) = log L(ω) L 0 = - j∈ω d j µ j . (1)
In real scenarii however, several sources of noise affect the image formation, and therefore a more comprehensive image formation model must be devised as illustrated in Fig. 1. We consider a sequence of images I = {I i } acquired from a volume discretised as a voxel grid with attenuations V = {µ j }. For each image I i we have:

1. A known projection and integration matrix P i composed of the coefficients d j obtained from motion capture. In the ideal case, we would have

P i V = I i .
We denote the projection matrix concatenation for all images P = {P i }. 2. A 2D image noise variance θ i accounting for the light source, the amplifier, and the imaging sensor. 3. Geometric noise, ie. the errors in the projection P i . This includes the inaccuracy in the motion and projection estimation as well as the deviation from purely rigid motion. It is modeled by a warping matrix F i .

Bayesian model

Our aim is to recover the 3D attenuation V given the absorbance image sequence I, ie. to invert the model described previously. For this purpose we rely on a MAP estimation to find the optimal solution in terms of attenuation and noise:

{V * , {F i } * , {θ i } * } = argmax V,{Fi},{θi} p(V, {F i }, {θ i })|{I i }), (2)
where, assuming statistical conditional independence between images given the attenuation model:

p(V, {F i }, {θ i })|{I i }) ∝ p(V ) i p(F i ) i p(θ i ) i p(I i |V, F i , θ i ). ( 3
)

Priors and image likelihood

Geometric noise appears as a result of calibration inaccuracies and non exactly rigid object motions. We modeled it by a warping function F i , estimated using the optical flow w i [START_REF] Bruhn | Lucas/Kanade meets Horn/Schunck: Combining local and global optic flow methods[END_REF] between the observed image I i and the generated one P i V .

As the inverse problem ( 3) is ill-posed and noise-ridden, we introduce noise and model priors. Given the nature of the data typically observed, the sparsity of the derivative responses is used as a prior for the 3D attenuation volume as in [START_REF] Liu | On Bayesian adaptive video super resolution[END_REF]:

p(V ) = η dim(V ) e -η ∇V , (4)
where η is the gradient weight. The minimisation of the L 1 norm of the gradient, or Total Variation T V L 1 , favours continuous volumes separated by potentially high, albeit localised gradients. The likelihood distribution is modeled as an exponential distribution:

p(I i |V, F i , θ i ) = θ dim(Ii) i e -θ i I i -F i P i V . (5)
where the 2D image noise variance θ i follows a Gamma distribution [START_REF] Liu | On Bayesian adaptive video super resolution[END_REF].

Model estimation

In order to solve for the parameters in the MAP estimation (2), we use a coordinate descent scheme [START_REF] Liu | On Bayesian adaptive video super resolution[END_REF], that iteratively cycles through the independent estimation of each parameter: the original volume V using expressions (4) and ( 5), image noise variance θ i , and the motion noise warp F i as detailed in section 3.3.

Attenuation Volume Given the current estimates for the warping function F i and noise θ i , we estimate the volume based on the image set {I i } and the gradient prior η with Iteratively Reweighted Least Squares (IRLS) to minimise: . Left-to-right: ground-truth CT scan, proposed method, without optical flow, without T V L1 prior, ART [START_REF] Pansiot | Combined visible and X-ray 3D imaging[END_REF]. Without T V L1 prior, the algorithm does not converge. The contrast is better with the proposed approach (better MI as compared to ART) even though artefacts appear on the edges as a result of aliasing during the data simulation process (higher RMS as compared to ART). ART performs relatively well in part due to the fact that simulated data are close to the noiseless theoretical model.

V * = argmin V i θ i F i P i V -I i + η ∇V . ( 6
Image/sensor noise Given the current volume V of N voxels and flow field F i we estimate the image noise variance θ i based on the residual error following a gamma distribution:

θ * i = α + N -1 β + x with x = N q=1 |(I i -F i P i V )(q)|, (7)
Geometric correction The residual motion is estimated using the optical flow w i [START_REF] Bruhn | Lucas/Kanade meets Horn/Schunck: Combining local and global optic flow methods[END_REF] between the observed image I i and the projected volume P i V . Given the current volume V and the noise variance θ i , we estimate the flow w i associated to the warp matrix F i . We then reformulate the data term in (6) as:

i θ i P i V -F -1 i I i . (8)

Experiments

Three sets of experiment were carried out to validate the proposed framework. First the CT scan of a phantom model was used to simulate image observations. Secondly, the phantom model was placed into our hardware platform. And third, an in-vivo hand was captured and reconstructed using the proposed framework.

Simulated radiographic and video data from CT

A forearm phantom, consisting of a real human forearm skeleton cast in resin was first scanned with a regular CT device. A complete capture pipeline was then simulated from the phantom scan which was rendered by 10 virtual video cameras and one virtual planar X-ray image 1 . The phantom scan model was then moved artificially, following roughly a 180 degrees rotation. The phantom motion was estimated using video and Iterative Closest Point for comparison with [START_REF] Pansiot | Combined visible and X-ray 3D imaging[END_REF]. The proposed approach was then applied to the simulated data. The performance of individual components of the algorithm were analysed independently, as illustrated in Fig. 2. We note that on simulated data, our approach exhibits slightly better contrast than ART (higher Mutual Information (MI) score). We also evaluate the sensitivity of our method with respect to the number of input frames, as illustrated in Fig. 3. These experiments show that for the given dataset, main skeletal structures can be recovered with as little as 16 frames. However, finer features such as bone cavities require at least 32 frames.

Furthermore we reduce the motion range from 180 to 90 degrees, which clearly impacts the volumetric estimation quality, as illustrated in Fig. 4. Raycasting rendering of the volume yields sharper results for poses within the original motion range, due to increased sampling density, as illustrated in Fig. 5.

In-situ forearm phantom

The proposed platform is composed of ten colour video cameras and a single X-ray C-arm. The forearm phantom presented here above was placed into the capture platform and moved manually to follow roughly a 180 degree rotation. of the phantom, as illustrated in Fig. 6. Furthermore, the energy spectrum of the CT scanner and that of the low-dose X-ray C-arm are different. Hence, the two models are first registered using multi-resolution Mutual Information (MI). The MI score is provided for quantitative comparison, being invariant to attenuation spectrum. Unlike the simulated case, this experiment shows that the proposed method performs significantly better than ART. In particular, the use of optical flow for motion noise compensation allows to retain a fair level of detail, whilst the T V L 1 norm prior constrains the ill-posed problem without excessive blurring.

In-vivo human hand

Finally, an actual in-vivo human hand was moved in the field, again following roughly a rotation movement over 20 frames. The results presented in Fig. 7 demonstrate the benefit of our approach which improves the results in some specific areas which we attribute to local (ie. non-rigid) motion. This demonstrates the interest of the generative model with the optical flow correction and the TV L1 regularization over more traditional approaches for few-view CBCT.

Conclusions and Future Work

In this paper we have presented a novel generative model to estimate a dense volumetric attenuation model of a rigidly moving object using motion tracking and a single planar X-ray device. Our framework contributes with an approach that takes benefit of object motion to accumulate evidence on the object inner structure. To this aim, we have introduced a Bayesian approach that optimally exploits X-ray information while enabling for acquisition noise. Our experiments show that the T V L 1 prior on the attenuation volume fundamentally contributed to convergence without excessive blurring, and that geometric noise can be effectively corrected using optical flow. This work considers rigid motion and we are currently investigating non-rigid motion.

Fig. 1 .

 1 Fig. 1. X-ray image formation model for a moving sample observed by a single static planar X-ray device. Video cameras are used to recover the sample motion.

Fig. 3 .

 3 Fig. 3. Results on simulated data (selected slice) based on varying numbers of input frames. Left-to-right: 8, 16, and 32 frames. Skeletal structures are visible with 16 frames when detailed features require 32 frames.

Fig. 4 .

 4 Fig. 4. Results on simulated data (selected slice) based on varying input angular range. Left: 32 frames roughly distributed over 180 degrees; right: over 90 degrees.

Fig. 5 .

 5 Fig. 5. Results on the simulated data (raycasting rendering) based on varying input angular range. Left: 32 frames roughly distributed over 180 degrees; right: over 90 degrees. The rendered viewpoint falls within the range of the original 90 degrees motion, but not on an original viewpoint, leading to sharper rendering due to locally denser sampling.

Fig. 6 .

 6 Fig.6. Results on the forearm phantom (2 selected slices and Mutual Information (MI) score). Left-to-right: ground-truth CT scan, proposed method, without optical flow, without T V L1 prior, ART[START_REF] Pansiot | Combined visible and X-ray 3D imaging[END_REF]. Without optical flow, artefacts are visible, for example in the bone cavities. The ART method produces much noisier results.

Fig. 7 .

 7 Fig.7. Results on the in-vivo hand with different priors (selected slice). Left-to-right: T V L1 weight η = 2; η = 1; η = 0 (no T V L1 prior); no optical flow prior; ART[START_REF] Pansiot | Combined visible and X-ray 3D imaging[END_REF]. The first three reconstructions demonstrate the favourable impact of the T V L1 prior. Comparing the first and fourth reconstruction (no optical flow), we observe less artefacts with flow and the bones are better resolved, in particular the index. ART exhibits a lot of under-sampling artefacts that cannot be recovered without priors.

  )

	RMS error	0.125	0.134	0.180	0.072
	MI score	0.338	0.318	0.177	0.309

Fig. 2. Results on simulated data (2 selected slices, RMS, and Mutual Information (MI) score)

Input data and more results for simulated and in-vivo experiments are available here: http://hal.inria.fr/hal-01348939/file/2016_Pansiot_MICCAI_Xrays3d_v2.mp4.

Acknowledgements

This research was partly funded by the KINOVIS project (ANR-11-EQPX-0024).


Baptiste Jonglez

Bruno Gaujal

Distributed Adaptive Routing in Communication Networks

Keywords: Multi-commodity routing, SDN, Nash equilibria, Adaptive routing Routage multi-flots, SDN, Équilibre de Nash, Routage adaptatif

In this report, we present a new adaptive multi-flow routing algorithm to select endto-end paths in packet-switched networks. This algorithm provides provable optimality guarantees in the following game theoretic sense: The network configuration converges to a configuration arbitrarily close to a pure Nash equilibrium. In this context, a Nash equilibrium is a configuration in which no flow can improve its end-to-end delay by changing its network path. This algorithm has several robustness properties making it suitable for real-life usage: it is robust to measurement errors, outdated information and clocks desynchronization. Furthermore, it is only based on local information and only takes local decisions, making it suitable for a distributed implementation. Our SDN-based proof-of-concept is built as an Openflow controller. We set up an emulation platform based on Mininet to test the behavior of our proof-of-concept implementation in several scenarios. Although real-world conditions do not conform exactly to the theoretical model, all experiments exhibit satisfying behavior, in accordance with the theoretical predictions.

Introduction

Communication networks are becoming increasingly multipath. Backbone networks have long been designed with multiple redundant paths, in order to increase reliability and performance. Datacenter networks need high bisection bandwidth for applications like MapReduce, which are typically obtained using multiple parallel paths in topologies such as fat-tree [START_REF] Al-Fares | A scalable, commodity data center network architecture[END_REF]. At the network edge, multihoming [START_REF] Abley | Goals for IPv6 Site-Multihoming Architectures[END_REF][START_REF] Troan | IPv6 Multihoming without Network Address Translation[END_REF] is increasingly common, due to the reduced cost of Internet access and prevalence of mobile networks such as LTE, as a complement to fixed landline Internet access.

A common challenge is to exploit this path diversity to increase performance and/or reliability. This problem can be modeled as a multi-commodity flow problem [START_REF] Hu | Multi-commodity network flows[END_REF]. Given a number of concurrent source-destination flows, the problem is to assign these flows to network paths, while respecting capacity constraints and optimizing a performance metric. Our goal is to provide a distributed solution to this problem, that requires neither cooperation between sites nor knowledge of the network topology and performance parameters. We restrict ourselves to the non-splittable case, in which a given source-destination flow uses only one path at a given time. This is necessary in practice to avoid persistent packet reordering, which would harm performance. However, a given flow is allowed to change its choice of path several times during its lifetime.

Related work

On the theoretical side, similar problems have been studied. When optimizing for a single flow with no a priori knowledge about the network, a multi-armed bandit modelisation [START_REF] Awerbuch | Adaptive routing with end-to-end feedback: Distributed learning and geometric approaches[END_REF][START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF] can be used for learning the shortest-delay path. However, since we consider multiple competing flows, we need a more complex modelisation based on game theory.

Regarding the multi-flow problem, convergence results for a wide class of online decision algorithms have been proven in [START_REF] Krichene | On the convergence of no-regret learning in selfish routing[END_REF], but in a different setting. The authors consider non-atomic routing games, i.e. each packet in a given flow independently chooses a path. This setting is hard to apply in a practical network because of packet reordering.

On the practical side, most existing approaches for exploiting multiple paths in a network are either limited to very specific settings, e.g. ECMP (Equal-Cost Multi-Path forwarding) for parallel paths with similar characteristics, or targeted at datacenter networks [START_REF] Cui | Difs: Distributed flow scheduling for adaptive routing in hierarchical data center networks[END_REF][START_REF] Wu | Dard: Distributed adaptive routing for datacenter networks[END_REF]. As such, they typically assume a network with a very organized and hierarchical topology, for instance fattree, or require explicit information sharing from the network [START_REF] Wu | Dard: Distributed adaptive routing for datacenter networks[END_REF]. We make no such assumptions.

When working with end-hosts directly, Multipath TCP [START_REF] Barre | MultiPath TCP: From theory to practice[END_REF] allows using multipath paths for a single source-destination flow, improving performance and reliability. This is the most promising approach so far, but requires extensive modifications in the operating system of both ends of the communication, which hampers large-scale deployment for now. In contrast, our proof-ofconcept only requires end-hosts to include TCP timestamps or use a LEDBAT-based transport protocol, which is already widespread [START_REF] Veal | New methods for passive estimation of tcp round-trip times[END_REF].

A noteworthy class of routing algorithms bases its routing decisions on dynamic properties of the network, including real time link load, end-to-end latency, or packet loss. This class of strategies is known as adaptive routing. Such solutions are attractive because they typically enable more efficient usage of network resources. In contrast, traditional load-oblivious routing algorithm may blindly over-utilize some links, leading to congestion, while other links have unused capacity. However, despite years of research, adaptive routing techniques did not see wide adoption. The main reason is the presence of potential instabilities and routing oscillations, which could make the cure worse than the disease. Indeed, early experiments with delay-based routing in the ARPANET resulted in severe stability issues under high load, rendering the network close to unusable [START_REF] Khanna | The revised arpanet routing metric[END_REF]. A possible tradeoff is to adapt routing paths at a much slower timescale. For instance, Link Weight Optimisation [START_REF] Fortz | Optimizing ospf/is-is weights in a changing world[END_REF] chooses paths based on a "representative" traffic matrix, which is typically updated no more than once a day. This solves the stability issue, but it is no longer possible to adapt to real-time network conditions. More recent solutions have been proposed for adaptive routing [START_REF] Cui | Difs: Distributed flow scheduling for adaptive routing in hierarchical data center networks[END_REF][START_REF] Kvalbein | Multipath load-adaptive routing: Putting the emphasis on robustness and simplicity[END_REF][START_REF] Wu | Dard: Distributed adaptive routing for datacenter networks[END_REF]. Some of these solutions are based on heuristics and offer no real guarantees on their performance. Furthermore, stability issues are typically overlooked, or no convergence guarantee is provided.

Contributions

In this work, we present a novel algorithm for adaptive routing in packet-switched networks, mapping source-destination flows to paths. We claim that it provides a viable and stable solution to adapt to traffic conditions, and effectively avoids congestion. Our algorithm is based on strong theoretical grounds from game theory, while our proof-of-concept leverages SDN protocols to ease implementation. Our routing algorithm is endowed with the following desirable properties for efficient implementation: -It is fully distributed: only local information is needed, and it requires no explicit coordination between routers; -It is oblivious to the network topology; -It is robust to outdated and noisy measurements; -There are no endless oscillations; -It does not require clock synchronization between routers, or between routers and end hosts.

We start by describing the high-level algorithm, assuming an idealized routing problem. This model assumes that the network allows source routing, i.e. end-hosts can choose to send packets along any network path to a destination. We then derive a practical algorithm that is equivalent to the first one, but readily implementable in a traditional next-hop forwarding paradigm. It operates in a more limited setting: a number of multihomed sites are connected at the edge of a core network. The algorithm then runs in the gateway of each site, which has multiple paths to the core network. In short, only gateways take routing decisions on behalf of flows, and the choice is limited to a choice of a next-hop router for each source-destination flow traversing the gateway.

We then describe our proof-of-concept implementation of this second algorithm, using Openflow to 1) program each router 2) collect feedback about the state of flows traversing the router. Note that our algorithm is fully distributed and does not require a central view of the network: as such, each router can be controlled by its individual SDN controller, typically collocated with the router itself.

Finally, we report and comment on experiments using Mininet [START_REF] Lantz | A network in a laptop: rapid prototyping for software-defined networks[END_REF] to evaluate the performance of our proof-of-concept in various scenarios.

2 Distributed Routing over a Network

Problem description and definitions

Let (V, E) be a communication network over a set V of nodes and a set E of bi-directional communication links, over which we consider the following multi-commodity flow problem.

A set K of flows of packets must be routed over the network. Each flow k ∈ K is characterized by a source-node a k , a destination-node b k and a nominal arrival rate of packets, λ k . Also, each flow is affected a set P k of possible paths in the network from its source to its destination, with Inria |P k | = P k . Notice that P k can either include all possible paths from a k to b k , or a single path, or any set in between these two extreme cases. A configuration is a choice of one path per flow.

Our objective is to find a configuration that minimizes a performance index. Here, we choose to minimize end-to-end delays of each flow. Actually, other criteria could have been chosen, such as loss rate, round trip times, or goodput, that would require a minimal adaptation of our algorithm.

This type of question has been heavily studied in the literature in many different ways [START_REF] Hu | Multi-commodity network flows[END_REF]. Our objective here differs from most previous work: We design a learning algorithm that allows each flow to discover a path, such that the global configuration is a Nash equilibrium of the system. Namely, after convergence, no flow can improve its delay by changing its path.

The challenge is to consider a realistic scenario in which no flow has information about the choices of the others or is even aware of the presence of other flows. The only information that a flow can get from the network is an estimation of the end-to-end delay of its packets sent over its current chosen path. To make things even more realistic, we assume that the delay measurements can be perturbed by random noise.

The difficulty also comes from the interdependence between the flows: when one flow decides to send its packets on a new path, this may alter the delays for the other flows because of resources sharing between the paths (either links or routers).

In this context we design a distributed algorithm that allows each flow k to eventually choose a path in P k such that all alternative paths would offer a larger delay. This goal is achieved by using optimization methods coming from game theory.

An illustrative example is given in Figure 1 with three flows. Flow 1 can be seen as the "main" flow while the other two flows are cross traffic that may alter the performance of the main flow.

We start from the configuration 1.(b), where all flows share some links, hence inducing large delays. The goal is to achieve one of configuration 1.(c) or 1.(d), by letting all flows explore simultaneously their paths and achieving such a coordination in a fully distributed way, in the sense that no flow has information about the presence of other flows and only has local information about the network.

OPS Algorithm

Here is a description of the algorithm that each flow k runs independently. It is based on a mirror-descent algorithm for general potential games, presented in [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF]. In the following, the index k may be dropped when no confusion is possible, but one must keep in mind that all flows execute the same code, not necessarily synchronously.

For one flow, say k, we call

d k (p 1 , • • • , p k , • • • , p K )
, the end-to-end average delay for packets of flow k under the configuration where flow 1 uses path p 1 among its possible paths, flow 2 uses path p 2 , and so forth.

The algorithm executed for flow k is probabilistic and maintains two vectors, both of size P k (denoted P to ease notations).

The probabilistic choice vector q = (q 1 . . . q P ) gives the probability to choose each path p.

The score vector Y = (Y 1 . . . Y P ) maintains a (negative) score for each path, to be optimized, where Y p depends on the average delay for packets of flow k on path p.

The main loop of the OPS algorithm (Algorithm 1) is as follows. Each time a local timer ticks (further details on this are provided in Section 3), a path p is chosen according to the probability distribution q, and packets are sent along p. The delay of packets over this path is measured. The score Y p is updated according to a discounted sum and in turn, the probability vector is modified for the next path selection using a logit distribution. This repeats forever, or

B 1 A 2 A 3 B 3 B 2 A 1 B 1 A 2 A 3 B 3 B 2 A 1 (a) Network. (b)Initial configuration. B 1 A 2 A 3 B 3 B 2 A 1 B 1 A 2 A 3 B 3 B 2 A 1 (c)Desired configuration.
(d)Desired configuration. until a stable path has been reached for all flows, i.e. q becomes a degenerate probability vector (all coordinates are equal to zero except one) for all flows.

Algorithm 1: Optimal Path Selection (OPS) Algorithm for flow k.

1 Initialize: 2 n ← 0; q ← ( 1 P , . . . ,

p ← Y p -γ n (D + τ Y p )/q p ∨ β n ; 9 foreach path s ∈ P k do 10 update probability: q s ← exp(Ys) exp(Y ) ;
11 until end of time;

The OPS algorithm uses 3 parameters: τ > 0 is a discounting factor over past scores. The bounding sequence β n (in the algorithm, ∨ denotes the maximum operator) is such that β n → -∞ and |β n | C 1 n+C 2 for some constants C 1 and C 2 . The decreasing sequence of discretization steps γ n is in L 2 ( n γ 2 n converges), but not in L 1 , ( n γ n diverges). Typically, γ n = 1/n α with 1/2 < α 1 works.

Discussion on how to tune these parameters is postponed to Section 3.6.

Inria

Convergence Properties

Assumption 1. The measurement D(t) of the delay is not biased. More precisely, the measurement done by flow k, over path p k at physical time t is such that D(t) = d k (p 1 (t) . . . p k (t) . . . p K (t))+ ξ(t), where the noise ξ(t) is a non-biased random variable, conditionally to the past, F t : E(ξ(t)|F t ) = 0, and has a finite second moment:

E(ξ(t) 2 |F t ) < ∞.
Recall that d k is the average end-to-end delay experienced by flow k in a given configuration. This is a rather mild assumption, because the errors on the measurements do not need to be bounded, nor does it have to be independent of the current configuration of the traffic, for example a large load over one link may induce a noise with large variance on that link. Assumption 2. All the flows have the same arrival rate.

This assumption is technical and deserves some comments. It is needed to make sure that the underlying game has a potential [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF]. This is an essential ingredient to prove convergence (see Appendix A).

However, it is not realistic: In general, all flows do not have the same arrival rate. In that case, one can still use the OPS algorithm. It can be shown (not reported here) that if OPS converges, it finds a Nash equilibrium. However, convergence is not guaranteed in general. An alternative, when all rates are commensurate, is to split flows into subflows, all with the same rate. In that case the subflows from the same flow may end up using different paths from source to destination. This solution is not ideal because it implies re-ordering of packets at the destination. However, if the arrival rate of the sub-flows is small enough, the delays on the different paths will be very similar after convergence, so that re-ordering will be minimal. The similarity of delays on different paths comes from the fact that, with small arrival rates, Nash equilibria approach Wardrop Equilibria, under which all the packets of a flow experience the same delay.

Theorem 1 (Convergence to equilibrium).

Under assumptions 1 and 2, for all > 0, there exists τ > 0 such that under discounting factor τ , Algorithm 1 converges for all flows to an -optimal configuration, in the following sense: For each flow k, the probability vector q converges almost surely to a near degenerate probability: q p becomes smaller than for all p ∈ P k except for one path, say p * k , for which it grows larger than 1 -. Furthermore, under configuration (p The proof of Theorem 1, given in Appendix A, also shows the following additional properties of the OPS Algorithm.

The algorithm is multiflow. The algorithm is distributed and multiflow in the sense that all the flows execute the algorithm in parallel and they all converge to a stable configuration.

The algorithm is incrementally deployable. Not all flows in the network need to use the OPS algorithm. If only a subset of the flows use the OPS algorithm, while all the other flows are forwarded using static routing tables, the convergence property still holds. The flows being forwarded statically are seen as part of the (random) environment over which optimization is done.

OPS is only based on online local information. The only information used by the flows is the current delay on their current path. They do not know the topology of the network, the capacity of the network, the number of concurrent flows or the paths taken by the other flows. The discovery of the best path for each flow is done "in the dark", without any exchange of information or any kind of collaboration between the flows or with the network.

Asynchronous updates. The algorithm is truly asynchronous: convergence will occur even if the different flows update their choices using different timers, that can have arbitrary initial offsets as well as arbitrary different speeds. Of course, the speed of convergence may be slow if flows have very different update frequencies.

Robustness to errors. As explicitly specified in Assumption 1, convergence to the equilibrium is robust to measurement errors. Actually, this is very useful in our context, not only because measurements do not have an infinite precision, but mainly because the measurement of the delay over one path is done by monitoring one or several packets (see Section 3.2 for a detailed description on how this can be done) and by computing the empirical mean. This empirical mean equals the average delay up to some additive random term. This term satisfies Assumption 1 as soon as the system is stationary.

Robustness to outdated measurements. The estimation of the delay D on path p is based on monitoring the delay of several packets taking this path. When the score Y p is computed and a new path is selected, these measurements are outdated because other flows may have already changed their paths. The theorem shows that this phenomenon does not jeopardize the convergence properties of the algorithm.

Convergence to a non-randomized choice. Theorem 1 says that there are no endless oscillations. The flows explore their different options in a transient phase of the algorithm. This creates some oscillations in the traffic. However, convergence to a stable point starts to settle in after some time and complete stabilization to the choice of a single path for each flow is always achieved.

Relaxing convergence

The guarantee of almost sure convergence requires that the step-size sequence γ n vanishes to 0.

One can relax this vanishing condition and show that for constant step sizes (γ n = γ, ∀n ∈ N), convergence occurs in probability instead of almost surely (see [START_REF] Benaïm | Dynamics of stochastic approximations[END_REF]). Namely, the distribution of q concentrates to near degenerate probability distributions when γ and τ go to 0.

Speed of convergence

The proof of the theorem given in Appendix A is based on the construction of a stochastic approximation of a differential equation and does not provide any information on the speed of convergence of the algorithm to a Nash equilibrium. The speed of convergence of similar algorithms has been proved to be of order 1/n [START_REF] Semen | Problem Complexity and Method Efficiency in Optimization[END_REF] in a centralized context (a single player) and with strictly convex objective functions. Neither conditions are true here (many players with arbitrary delay functions). To our knowledge, no theoretical result exists today to bound the speed of convergence for the type of algorithms used here, so it must be evaluated experimentally.

From a practical point of view, several factors can affect the convergence of the OPS algorithm.

• Size of the network: The algorithm being based on the measurement of end-to-end delays, the size of the network should not play a decisive role in the convergence speed. The effects of the network size are indirect: a larger network means a larger number of possible paths for each flow, and more potential for delay measurements to be outdated.

• Number of players: Classical algorithms let players play one at a time, while our OPS algorithm makes them play simultaneously. This speeds up convergence, especially with a large number of players. Several simulations of OPS done in a different context (not reported here) suggest that the speed of convergence does not depend crucially on the number of players.

Inria

• Number of possible paths per player: In contrast with the number of players, the convergence time is expected to increase drastically when the number of choices per flow increases. Each flow needs to probe each path a sufficient number of times before it can assess its performance.

• Precision of the measurements: It should be clear that when the measurements of the endto-end delays suffer from a large variance, then convergence slows down. We believe that this is the most important parameter that influences the convergence time. This belief is reinforced by the experiments reported in Figure 5.

Price of Anarchy

In general, a Nash equilibrium (NE) does not provide any guarantee on its global performance.

In the worst case, the sum of the delays of all flows under a Nash equilibrium can be arbitrarily far from an optimal configuration. However, there are more favorable situations in which Nash equilibria can provide performance guarantees.

Here, the social cost of any path configuration p = (p 1 , . . . , p K ) is the total sum of the average delays that can be decomposed over all links of the network, c(p

) := e∈E N e (p) • δ e (N e (p)),
where δ e ( ) is the average delay on link e under load and N e (p) is the number of flows using link e under configuration p.

A first argument in favor of NE is that social optima that are not NE are unstable configurations, because some flows can gain by changing their paths, and will certainly do so if they are not constrained by some kind of central authority.

Moreover, in some cases, NE exhibit good global performance. This has been deeply investigated in the literature: we simply provide a quick review here. In general, the price of anarchy (ratio between the social cost of the worse Nash equilibrium over the cost of the social optimum) is upper-bounded by Θ(log n) [START_REF] Christodoulou | The price of anarchy of finite congestion games[END_REF] where n is the number of players.

Furthermore, if the delays of all links are (λ, µ)-smooth (see the definition in [START_REF] Christodoulou | The price of anarchy of finite congestion games[END_REF]), then the price of anarchy is bounded by a constant, namely λ 1-µ . As a special case, this bound becomes 5/2 when the delay on each link of the network is affine w.r.t. its load (see [START_REF] Christodoulou | The price of anarchy of finite congestion games[END_REF]). Additionally, when the delays are arbitrary increasing functions of the load, the social cost of any NE is smaller than the optimal social cost of the same congestion game where rates are doubled [START_REF] William | Population Games and Evolutionary Dynamics[END_REF].

There exists cases where NE are local social optima. This is for example the case when the potential function (minimized by NE), φ(p) = e∈E Ne(p) i=0 δ e (i) approaches the social cost c(p). This can happen when the delays do not depend on the load of links (M/M/∞ models for the delay on all links). However, this may not be the most interesting case, because then the flows will not interfere with each other at all.

Finally, several recent publications question the worst case bounds found in [START_REF] Roughgarden | How bad is selfish routing[END_REF] and show that in real life cases with a large number of players [START_REF] Cole | The price of anarchy of large walrasian auctions[END_REF] or when the arrival rate grows [START_REF] Colini-Baldeschi | On the price of anarchy of highly congested nonatomic network games[END_REF], the price of anarchy approaches one.

Implementation as a SDN controller

We now focus on the more practical aspects of our model, and discuss how Algorithm 1 can be implemented in a real network. We then describe a more practical version of our algorithm, designed as a SDN controller. We have successfully implemented this SDN-based algorithm, and we use this implementation to run experiments in Section 4.

Running our algorithm in gateways

Algorithm 1 assumes that flows, or equivalently end-hosts, can choose the full path to their destination (line 6 of the algorithm). This paradigm is known as source routing, but is not widely deployed in the public Internet: source routing was deprecated in IP because of security issues [RFC 5095]. Thus, we first make assumptions on the structure of the network, and we then adapt Algorithm 1 so that routing decisions are taken by routers on behalf of end-hosts.

First, we assume that the network can be decomposed into a core network and several edge networks. Each edge network connects to the rest of the network through a single router, which we call a gateway. Gateways may connect directly to each other, or through core routers. Figure 2 shows an example with three edge networks, each of which contains a gateway and one or more hosts.

Then, to ease implementation, we assume our OPS algorithm is only run by gateways to select a path among those offered by the core network. This way, core routers do not need to keep state for each source-destination flow. Core routers simply use their regular routing protocols, for instance BGP or OSPF. They can also use ECMP (Equal-Cost Multiple Paths) with a hashing mechanism on the packet source and destination, so that a given source-destination flow in the core network is forwarded along a unique path. This is necessary to obtain consistent delay measurements during the lifetime of a flow.

Limiting our algorithm to only run in gateways means that a flow cannot explore all possible paths to its destination: the set of possible paths P k defined in Section 2 for flow k is limited to the next-hop routers of the gateway. Still, even with two or three paths to choose from, this provides enough path diversity to make our algorithm interesting to run.

The decomposition of the network into gateways and core routers is illustrated in Figure 2.

Exploiting one-way delay

Algorithm 1 uses a generic notion of "delay" as objective function (line 7 of Algorithm 1). We focus on end-to-end one-way delay instead of the more classical RTT. There are several reasons to this choice: Inria 1) For latency-critical applications such as Voice-over-IP, the relevant metric is often the transmission delay on the forward path. 2) A router has no control over the reverse path since routing decisions only affect the forward path of a flow towards its destination. For instance, the LEDBAT congestion control algorithm [START_REF] Shalunov | Low Extra Delay Background Transport (LEDBAT)[END_REF] also exploits one-way delay for the same reason: it allows congestion detection and bufferbloat avoidance on the forward path, while ignoring the effect of the reverse path.

To passively determine the one-way delay of forwarded packets, we use timestamps that some transport protocol include in their headers: we focus on TCP with the Timestamp Option and µTP, an implementation of LEDBAT [START_REF] Shalunov | Low Extra Delay Background Transport (LEDBAT)[END_REF] over UDP. To estimate one-way delay from these timestamps, we use the same method as [START_REF] Chirichella | Passive bufferbloat measurement exploiting transport layer information[END_REF] to extract one-way delay samples from TCP and µTP/UDP packets, and then normalize the measurements by expressing them in milliseconds.

In our case, the gateway is the "observation point" at which the delay is observed. A gateway is interested in the one-way delay from A to B, where A is a host in the local edge network, and B is a remote host. The observed one-way delay is measured for every packet i from A to B, and we call it δ observed (i). It is the only quantity we can measure directly, but it can be theoretically decomposed in this fashion:

δ observed (i) = δ propagation (p i ) + δ queuing (i) + off A,B
where p i is the path taken by packet i in the network. Here, we consider three contributions to the observed one-way delay:

1. the clock offset of B relatively to A, off A,B . This quantity is independent of the path taken by a packet, and we further assume that it does not change over time.

2. a delay contribution that depends only on the path taken by a packet, δ propagation (p i ). This is the one-way propagation delay defined by physical laws (e.g. the speed of light in a fiber) and by the constant switching time of switches and routers on the path. By definition, this quantity is the same for all packets forwarded on the same path.

3. the remaining contribution, which is specific to each packet. We call this last contribution δ queuing (i), because it is mostly made of queuing delays at switch and router interfaces along the path. This contribution actually also includes serialization delay, although it is negligible in high speed networks.

To eliminate the clock offset and get positive estimates of the delays (needed to guarantee proper convergence of the algorithm), we compute the minimum of the observed delay measurements, similarly to LEDBAT. In our case, this minimum is taken over all paths selected by the algorithm. Given a measurement of the observed delay δ observed (i) on the ith packet, with p i being the path selected for this packet, the relative delay we consider for our algorithm is the following:

∆(i) := δ observed (i) -min 0 l i δ observed (l) (1)
= δ observed (i) -min

0 l i δ propagation (p l ) + δ queuing (l) + off A,B (2)
= δ observed (i) -off A,B + min path p δ propagation (p) (3) = δ propagation (p i ) + δ queuing (i) -min path p δ propagation (p). (4)
Going from (2) to (3) requires that i is large enough to achieve a minimum queuing delay of zero, and to ensure all paths have been visited. In practice, our algorithm starts with a uniform probability to select each path, so most paths are expected to be probed before the algorithm starts to converge. This way, we eliminate the clock offset and obtain a positive estimate of the one-way-delay on the current path p i , relatively to the path with the lowest propagation delay. It can be shown that this additive shift does not affect the behavior of the algorithm, given that ∆(i) remains non-negative.

Update policy

Algorithm 1 (line 4) uses a timer for its updates, denoted n. This becomes a periodic update in the implementation: for each flow forwarded by a gateway, the gateway takes a routing decision every T packets of this flow. Among these T packets, the gateway only uses the one-way delay of the last S of them and compute their empirical mean. Using the last packets of the period helps to ensure that the transient state associated with the last change of path has disappeared, and that measured delays reflect the steady-state performance of the flow.

Recall that ∆(i) is the relative delay of packet i, defined in Section 3.2. Using the above notations, the delay D used in line 7 of Algorithm 1 is computed as:

D := ∆(nT -S + 1) + • • • + ∆(nT ) S .
This choice has several advantages:

• Large flows (i..e. with large arrival rates) are more likely to affect the delays of other flows. Using this construction of D makes sure that they are updated more frequently.

• This construction also provides some control on the random error ξ (defined in Assumption 1): A large T ensures that the expectation of ξ is close to 0 while a large S reduces its variance.

SDN-based implementation

We are now ready to adapt Algorithm 1 so that it is more practical for real networks. Below is a description of the main ideas of our implementation, while Section 3.5 contains a more formal description of the practical algorithm we implemented. The code of our actual implementation is available online1 . Our proof-of-concept implementation takes the form of an Openflow controller, using the Ryu [START_REF]Ryu: a component-based software defined networking framework[END_REF] library. A gateway is made of both an Openflow switch and a dedicated controller, but we use the term "gateway" to refer to the Openflow switch only.

The forwarding table of a gateway is programmed and constantly updated by its controller. Furthermore, a gateway sends a copy of packet headers back to its controller, for delay computation. The controller is configured beforehand with a static multipath routing table. That is, for each IP prefix, the routing table lists all possible next-hop routers that can be used to reach the destination.

When a gateway receives the first packet of a new flow, it forwards it to the controller, to receive instructions for subsequent packets of the same flow. The controller takes an initial routing decision for this flow, according to its routing table. This decision is a choice of a nexthop router through which packets will be forwarded. A corresponding forwarding rule is then installed into the gateway.

To measure the performance of each decision, the controller also installs a rule to receive a copy of all packet headers. Upon receiving a header, the controller extracts timestamps, computes the end-to-end one-way delay of the flow, and updates the score of the current choice. For each active flow independently, the controller periodically decides to select a new next-hop router, based on the scores. Each time the controller changes its decision, new forwarding rules are installed in the gateway.

The main reason to use a SDN framework is ease of implementation. Indeed, Openflow abstracts away the communication with forwarding elements, allowing to easily install forwarding rules and receive packet headers for performance measurements. For our experiments, we used Openvswitch on Linux, but any Openflow-compatible switch could have been used instead.

An important architectural decision was to exploit the decentralized nature of our algorithm, by ignoring the centralized management features of Openflow. Here, each gateway is controlled by a separate Openflow controller. Thus, we retain the scalability and resilience properties of decentralized routing.

Another key point is that we use Openflow switches as pure layer-3 devices. Openflow 1.3 is perfectly capable of programming a switch to act as a layer-3 router, by decreasing the TTL and modifying the MAC addresses of packets. Working at layer 2 would not make much sense, because we fundamentally solve a routing problem. Besides, it is difficult to have control on flooding when working at layer 2 in complex topologies with redundant paths.

Practical algorithm

Algorithm 2 is a version of Algorithm 1 that is much closer to the actual implementation. This algorithm is run independently for each gateway, where a flow is identified by its source and destination addresses. That is, the controller of each gateway executes this algorithm independently, for all flows currently forwarded by the gateway.

Parameters

Both Algorithm 1 and Algorithm 2 use a number of parameters. We discuss how to choose appropriate values for these parameters, and the default value we use in our implementation. These default values are then used in the experiments we present in Section 4.

Step size of score updates γ (Section 2.2) We choose a constant value of γ, instead of a decreasing sequence γ n . This implies convergence only in a weaker sense (in probability instead of almost surely), as discussed at the end of Section 2.3. However, a constant value of γ is interesting in practice, because it reduces convergence time: in the original algorithm, when γ n goes to zero, convergence to the equilibrium becomes slower and slower as n increases. In short, we trade a bit of theoretical guarantees for a better convergence time.

The numerical value of γ should be small compared to typical measured delays, since it is the discretization step of the differential equation driving the dynamics of the system (see Appendix A). In our implementation, we have chosen γ = 0.01, which satisfies the above constraint since we use delays expressed in milliseconds.

Discount parameter τ (Section 2.2) A smaller value means that the algorithm will converge closer to a Nash Equilibrium, but convergence will be slower. We choose τ = 0.1. This allows fairly quick convergence, at the expense of a little bit of accuracy. Update interval T (Section 3.3) By default, gateways update their choice every T = 500 packets for each flow. Assuming a flow with a throughput in the range of a few Mbit/s, this means our algorithm performs an update every few seconds at most. Note that an update of the algorithm does not necessarily mean that a new path is selected: the new choice can be identical to the previous choice. Thus, the rate of actual routing changes will be much lower than this, especially since the algorithm is expected to converge to a deterministic choice.

A few seconds between choices is a reasonable value: the network has time to adapt before measuring the performance of the new configuration. Indeed, just after the change, the queues of routers in the network might exhibit a transient behavior, which is not adequate to measure average delays. Additionally, we measure delays using packets on the reverse path. If the endto-end RTT of a flow was higher than our update interval, then each decision would be based on outdated information. The algorithm is designed to be robust to such outdated information, but we still want to avoid it, as it may slow down convergence.

For flows with a larger throughput, T = 500 packets would be too low, as it would mean very frequent changes of path. A possibility would be to specify an update interval in seconds, independently of the throughput of the flow. However, this might cause fairness issues between flows with different throughput. Since we explicitly assume that all flows have the same throughput to be able to prove convergence, we have not yet explored this area.

Summary of parameters

The various parameters used in the implementation and their default value are summarized here: γ τ T S 0.01 0.1 500 packets 5 packets

Experimental Evaluation

We next evaluate our OPS algorithm in a simple network, to make sure it discovers optimal paths in practice. We are also interested in the time it takes to reach a solution.

The experiments were performed using emulation on Linux systems. This means that we were able to run our SDN implementation unmodified: we thus expect our results to be quite close to reality. Using emulation also means that we were unable to perform a large number of experiments or use large-scale topologies. Indeed, setting up and then running each experiment in a reproducible manner requires writing the configuration of each router and controller. We used Mininet [START_REF] Lantz | A network in a laptop: rapid prototyping for software-defined networks[END_REF] to emulate network topologies, running on a single server. Although realistic in a functional sense, Mininet is known to exhibit unrealistic performance characteristics [START_REF] Lantz | A network in a laptop: rapid prototyping for software-defined networks[END_REF]. To prevent this effect, we limit the capacity of the links to very low values using netem (in the order of 10 Mbit/s). Furthermore, we run our Mininet network on fast and modern server hardware, which should alleviate performance concerns given the low throughput involved in the experiments.

Experimental setup

To evaluate our algorithm, we consider a network satisfying the conditions laid out in Section 3.1. This network is simple enough to understand how the algorithm behaves, but complex enough so that the solution found by the algorithm is non-trivial. The network is shown in Figure 3. A set of gateways G 1 to G 4 are connected to each other, and to a central router R that models the core of the network. A set of hosts, each attached to a gateway, send long-lived UDP flows to each other. Five UDP flows are considered, numbered from 1 to 5: from Host 1 to Host 3, from Host 4 to Host 1, from Host 4 to Host 2, from Host 4 to Host 3, and finally from Host 2 to Host 3.

For this evaluation, all flows have a constant throughput equal to λ, and have an infinite duration. It can be seen as a model of a particular application usage, for instance a network used to transport multiple real-time video streams.

Note that our OPS algorithm only runs in the four gateways. The central router R simply forwards packets using a static routing table, designed to minimize the number of hops to the destination. In this case, all gateways are one hop away from R, so the routing table trivially forwards packet directed to host i through gateway G i .

Additionally, each gateway uses static routes for packets coming from other gateways since the OPS algorithm is only run for packets from its own hosts (see Section 3.1). The next-hop router of these static routes is the destination gateway if it is a neighbor, and R otherwise.

The capacity of each link is detailed in the Figure . To ease explanations, we introduce the load ρ as ρ = λ 4000 Kbit/s . A load of 1 means that a single flow would occupy the full capacity of a 4 Mbit/s link.

As long as ρ < 1, there exists at least one stable configuration of the network, i.e. a choice of path for each flow that satisfies capacity constraints on all links. Additionally, for 2 3 < ρ < 1, there are only two possible stable configurations, described in Figure 4.

Since we have four gateways, we run four instances of our Ryu-based Openflow controller, each controlling a different emulated gateway. Once Mininet and the controllers are setup, we then use udpmt (part of the ipmt toolbox [START_REF] Berger | ipmt: Internet protocols measurement tools[END_REF]) to generate UDP packets for each flow. We modified udpmt so that it sends packets with a µTP header, which includes a timestamp of the date of emission. The destination host for each flow runs udptarget, which we modified to behave similarly to a µTP implementation: it replies back with small UDP packets containing µTP timestamps.

This setup has been run a large number of times, to reduce variability. Each experiment lasts for 2400 seconds, because we found it was sufficient to allow most executions to converge. To parallelize the execution, we used Grid'5000 [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF] as an IaaS platform. We reserved several identical physical machines from a cluster, deployed a Debian Jessie image with all necessary software, and ran the same experiments independently on each machine. The machines are Bullx Blade B500 servers, with a 8-core Intel Xeon E5520 CPU and 24 GiB of RAM. Depending on the experiments, we used between 5 and 45 machines in parallel.

Flow Equilibrium 1 Equilibrium 2 Flow 1 G 1 → G 2 → G 3 G 1 → G 2 → G 3 Flow 2 G 4 → G 1 G 4 → R → G 1 Flow 3 G 4 → R → G 2 G 4 → R → G 2 Flow 4 G 4 → R → G 3 G 4 → G 1 → R → G 3 Flow 5 G 2 → G 3 G 2 → G 3

Discussion of the results

The experiment was carried under several values of load ρ (ρ = λ 4 Mbit/s ). In the experiments, the load varies from 0.75 to 0.975, to exercise the algorithm under moderate and heavy loads.

Figure 5 shows the evolution of Flow 1 under different loads. Gateway G 1 has three possible next-hop routers: G 2 , R, or G 4 . We display the probability, over time, that G 1 selects each of the three next-hop routers. The other flows are also being forwarded concurrently, but this is not displayed on the figure .  G 2 gets selected most of the time, after a transient period. This choice is consistent with both optimal configurations for the network. We also note that, when the load increases, the transient period becomes longer. Looking at Figure 5c, the convergence of Flow 1 is not obvious when the load is high. In fact, Flow 1 has converged to an "almost pure choice" as predicted by Theorem 1. This means that the flow will spend most of its time on the optimal choice, but can explore other choices from time to time. This effect is more visible when the load is high: in Figure 5c, G 2 is selected most of the time, but the flow also explores G 4 and R from time to time. Taking RR n°8959 q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 0 200 400 600 800

Number of iterations

q q q R G2 G4
Selection probability for each next-hop router (a) Load equal to 0.75 q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

q q q R G2 G4
Selection probability for each next-hop router (b) Load equal to 0.875 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 Inria into account these observations, we have designed the following convergence criterion for a flow:

A flow has converged to path p at time t if the average probability of choosing p for the last W updates was at least Θ.

Here, W is the size of a moving window counted as a number of updates of our algorithm (recall from Section 3.3 that a gateway periodically updates the path of its flows, every T packets).

The global convergence criterion for a network with multiple flows is satisfied when all flows satisfy the convergence criterion at time t, for the same parameters W and Θ.

Using this convergence criterion over the previous experiments with W = 50 and Θ = 80%, we have computed the global convergence time of the algorithm under different loads. Note that the convergence time cannot be lower than W = 50 iterations. The results are depicted in Figure 6, where the convergence time is expressed as a number of iterations of our algorithm. For a moderate load, the global convergence time is close to the minimum of W = 50 iterations. When the load approaches 1, we expect that the global convergence time goes to infinity, because the variance of the delays goes to infinity.

Note that we intentionally exercise our algorithm in unfavorable conditions. For a low or moderate load, which is the usual case for over-provisioned operator networks, the global convergence time is very close to the minimal convergence time.

Quality of the final configurations

In the network of Figure 3, there exist two stable configurations (in other words Nash equilibria) where no flow can lower its delay by changing its path (these two configurations are given in Figure 4). By Theorem 1, the OPS Algorithm converges to either of these configurations. Only Flows 2 and 4 differ between these two configurations. In the second configuration, both flows choose a longer path. Thus, the first configuration is better than the second one, because the average delay of any flow is either smaller or unchanged. However, both configurations successfully avoid congestion so that they both constitute good choices (even if not optimal).

The following table shows the percentage of experiments that converged to equilibrium 1, equilibrium 2, any non-stable configuration, or did not converge. As the load increases, it can be observed that equilibrium 1 (the best configuration) is more likely to be selected by our algorithm, which matches intuition. For each value of network load, 85 independent experiments were conducted. The results in the table use our convergence criterion with Θ = 0.8 and W = 50. The executions reported as non-convergent had failed to converge after the end of the experiment (2400 seconds).

The quality of the convergence actually depends on the parameters Θ and W . If one chooses the threshold Θ close to one, then the quality improves, in the sense that the number of times that the OPS algorithm stops under an unstable configuration drops. However this comes at a cost: the time of convergence increases as well as the number of executions that did not converge. An empirical search (not reported here) with varying values of Θ and W shows that a good compromise in our experiments is obtained by choosing Θ = 0.8 and W = 50.

Conclusion

We have described an adaptive routing algorithm for packet-switched networks, first as a theoretical algorithm, and then as a practical implementation that can be used as a Openflow controller. We have first proved a convergence theorem in a idealized network model, showing that our algorithm converges to a quasi-Nash Equilibrium. Using our prototype implementation in an emulated network, we have then shown that our SDN implementation converges fast under light to moderate network load, and indeed reaches a Nash Equilibrium (i.e. a stable configuration where no flow can improve its delay). Under very heavy network load, convergence is slower, which is not surprising given the variability of delays when the network becomes completely saturated. However, even in these extreme conditions, our implementation still converges to a quasi-Nash Equilibrium in most cases.

6 Future work

Dynamic network environment

Our current implementation uses static routing tables. An obvious improvement would be to couple our controller with a routing protocol, such as BGP or OSPF, that would compute possible paths for each destination and feed them to our algorithm.

This brings up the question of dynamic network conditions, which would create new paths or invalidate existing paths in the network. In addition, flows can have varying demand over time. For now, we have focused on routing over a stationary network. If the network becomes dynamic, the proof of convergence would be more difficult, and can no longer be based on the theoretical tools we used. One promising direction is to use two time-scale stochastic approximations: On the slow time scale, network conditions and flow demand can evolve. On the fast time-scale, the algorithm updates its scores and path choices for each packet.

Congestion control

Let us consider a more challenging situation, where flow endpoints use a congestion control algorithm. This means that the throughput may depend on the routing choices, since the packet loss, available capacity and delay for a flow all depend on the path used by the flow.

This raises the question of the interaction between congestion control at endpoints and routing decision made by our algorithm. This can possibly lead to a "self-confirming" situation as described by [START_REF] Eric | On the stability of adaptive routing in the presence of congestion control[END_REF], where congestion control reduces throughput because of bottlenecks created by poor routing choices, leading to this stable but inefficient situation. In [START_REF] Eric | On the stability of adaptive routing in the presence of congestion control[END_REF] the authors further prove that an adaptive routing algorithm that minimizes the maximum link utilization (the so-called "min-max" property) converges to a stable and optimal situation, given some fairness assumptions on the congestion control mechanism.

Thus, we leave to future work the stability and optimality considerations related to congestion control; for instance, one could study situations and hypotheses for which our algorithm satisfies the min-max property.

Figure 1 :

 1 Figure 1: (a) displays a network with 3 flows (A 1 , B 1 ), (A 2 , B 2 ), (A 3 , B 3 ). (b) shows the initial configuration that suffers from congestion in two links (marked by star shapes). The final configurations (c) and (d) have no congestion.

Figure 2 :

 2 Figure 2: Illustration of the network decomposition into hosts, gateways and routers. An example of 3 flows between hosts is also displayed.

5 if pkt belongs to an unknown flow then 6 flow1 7 take_routing_decision (flow1) 8 flow2 9 take_routing_decision (flow2) 10 else 11 Find 16 D

 56789101116 ← create_flow (pkt.src, pkt.dest) ← create_flow (pkt.dest, pkt.src) flow f such that f.src = pkt.dest and f.dest = pkt.src 12 ∆ ← compute_onewaydelay (f, pkt) 13 Store ∆ in a queue for later use 14 f.packet_count += 1 15 if f.packet_count T then ← average of the last S values of delay ∆ for flow f 17 update_score (f.current_nexthop, D) 18 take_routing_decision (f) 19 f.packet_count ← 0 // Parameters S and T are presented in Section 3.3 // Functions in boldface are detailed in Algorithm 3Bounding sequence β n (Section 2.2) This sequence is essentially useful for the proof and does not play an important role in practice. It is skipped in the implementation.

1 Gateway G 4 Gateway G 3 Gateway G 2

 1432

Figure 3 :

 3 Figure 3: Network used in the following experiments. Thin lines represent links with capacity of 4 Mbit/s, while thick lines represent links with a capacity of 8 Mbit/s. All links have a transmission delay equal to 5 ms. Dashed lines between hosts and their gateway are links with unrestricted capacity. Five flows using the network are represented. Each flow consists in UDP packets with a constant throughput λ, with λ varying from 2000 to 3900 Kbit/s in the experiments.

Figure 4 :4000 Kbit/s satisfies 2 3 ρ < 1 .

 431 Figure 4: Stable configurations for the network of Figure 3, when the load ρ =

  each next-hop router(c) Load equal to 0.925

Figure 5 :

 5 Figure 5: Probability of selecting each next-hop router over time, for Flow 1 (see Figure3). Each figure corresponds to an increasing value of load for all flows. Here, T = 500 packets and S = 5 packets.

Figure 6 :

 6 Figure 6: Average global convergence time of the algorithm as a function of the load, with 95 % confidence intervals.

  Algorithm 2: OPS Algorithm for the SDN controller Install a table-miss Openflow rule in switch: priority=0, match=*, action=CONTROLLER 4 When packet pkt arrives at the controller

	1 Initialisation
	2	Populate routing table with static entries (mapping IP prefixes to lists of next-hop
		routers)

3

RR n°8959

https://gforge.inria.fr/projects/derouted Inria

Inria

A Appendix: Proof of Theorem 1

This appendix is devoted to the proof of Theorem 1. The OPS algorithm is a variant of an algorithm presented in [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF], that computes Nash equilibria in potential games. First, our routing problem can be seen as a game. The flows are the players, their strategies are the choices of paths and the payoffs (or the costs, here) are the expected delays on the chosen paths. Under this framework, our problem fits in a well-known class of games, namely atomic routing games. In [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF], it was shown that atomic routing games where each player (flow) has rate 1 are potential games, i.e. there exists a potential function for such games, namely φ(p) = e Ne(p) i=1 δ e (i).

Here, the rates are not necessarily equal to one but they can easily be replaced by flows with rate one by changing units. This is where Assumption 2 is used in our case. It implies that our problem is a potential game.

Using this construction, our algorithm is similar to Algorithm 1b in [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF] with the addition of bounding terms β n .

The L2-L1 condition on γ n is the same as assumption A1 in [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF]. Our Assumption 1 implies Assumption A2 in [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF]. As for assumption A3 in [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF], it is replaced here by the explicit bounds β n . The bounded approximation theorem (Theorem 2) in [START_REF] Andradóttir | A stochastic approximation algorithm with varying bounds[END_REF] allows us to state that the sequence Y n is an asymptotic pseudo-trajectory (APT) in the sense of [START_REF] Benaïm | Dynamics of stochastic approximations[END_REF] as soon as the clock ticks for the flows all have finite rates. This implies that for each flow the scores of paths, Y p , approach the solution of the following differential system ( 5)-( 6), even if the estimation D is based on outdated measurements of packet delays

as long as this solution is locally stable (in the dynamical system sense).

In [START_REF] Coucheney | Penalty-Regulated Dynamics and Robust Learning Procedures in Games[END_REF], is is shown that the solutions of this differential system are locally stable when the game is a potential game and are (τ ) approximations of Nash equilibria of the game, where the error (τ ) vanishes as τ → 0. This means that Theorem 2 in [START_REF] Andradóttir | A stochastic approximation algorithm with varying bounds[END_REF] can be used here.

Finally, the existence of a potential function further implies that the only locally stable Nash equilibria are pure (see [START_REF] Monderer | Potential games. Games and economic behavior[END_REF][START_REF] William | Population Games and Evolutionary Dynamics[END_REF]). In our case this means that there exist a path p k for each flow k such that d k (p 1 , . . . , p k , . . . , p K ) d k (p 1 , . . . , p , . . . , p K ), for all p in P k . Therefore our algorithm will converge arbitrarily close to such a pure Nash equilibrium. This concludes the proof.

Inria