Jean-Michel Fourneau

Bruno Gaujal

Marc Lelarge

Jean Mairesse

Laurent Truffet

Bruno Tuffin

Sara Alouf

Alain Jean-Marie

Valeria Neglia

Laurie Vermeersch

Soutiens

Stéphane Durand
email: stephane.durand@ens-lyon.fr

M Estécahandy
email: maider.estecahandy@total.com

L Bordes

S Collas

C Paroissin

Marion Dalle
email: marion.dalle@cea.fr

Jean-Marc Vincent

Florence Perronnin

Damien Dosimont

Guillaume Huard

Hadrien Hours

Ernst Biersack

Patrick Loiseau

Felipe Olmos

Bruno Kauffmann

Alain Simonian

Yannick Carlinet

Nicaise E Choungmo Fofack
email: nicaise.choungmofofack@orange.com

Seregina Tatiana
email: seregina@laas.fr

M Larra Ñaga

U Ayesta

I M Verloop

Luka Stanisic

Arnaud Legrand

Farah Ait Salaht
email: farah.ait-salaht@prism.uvsq.fr

Konstantin Avrachenkov
email: k.avratchenkov@inria.fr

Natalia M Markovich
email: markovic@ipu.rssi.ru

Jithin K Sreedharan
email: jithin.sreedharan@inria.fr

A perfect sampling algorithm of random walks with forbidden arcs

Objectif

L'Atelier en Évaluation de Performances est une réunion destinée à faire s'exprimer et se rencontrer les jeunes chercheurs (doctorants et nouveaux docteurs) dans le domaine de la Modélisation et de l'Évaluation de Performances, une discipline consacrée à l'étude et l'optimisation de systèmes dynamiques stochastiques et/ou temporisés apparaissant en Informatique, Télécommunications, Productique et Robotique entre autres.

La présentation informelle de travaux, même en cours, y est encouragée afin de renforcer les interactions entre jeunes chercheurs et préparer des soumissions de nouveaux projets scientifiques. Des exposés de synthèse sur des domaines de recherche d'actualité, donnés par des chercheurs confirmés du domaine renforcent la partie formation de l'atelier.

Historique

Démarré sous l'impulsion de l'équipe de Raymond Marie en 1991, l'atelier en évaluation de performance a eu lieu à Rennes (1991-93), Grenoble (95), Versailles (96), Paris ENS (01), Reims (03) et Aussois (08). L'objectif est de fédérer une communauté allant des probabilités aux expérimentations en informatique. Elle recouvre donc d'une part les thématiques portant sur la modélisation des systèmes complexes avec les méthodes théoriques récentes, les environnements logiciels d'évaluation de performances jusqu'aux retours d'expérimentations en vraie grandeur.

Illustrations with Call Centers Résumé : Large systems that involve humans (such as health-care systems, call centers, emergency systems, transportation networks, supply chains, communication systems, etc.) are difficult to manage because they are complex and involve significant uncertainty which is itself hard to model in a realistic way. For example, call arrivals in call centers follow stochastic processes whose rates are themselves random and depends significantly on the time of the day, type of day (day of the week, holiday), period of the year, weather, other external events, etc. The arrival processes of different call types may also be dependent. Call durations (service times) have distributions that depend on the call type and on the particular agent who handles the call, and are often time-dependent because the effectiveness of agents depends on their experience, base qualities, motivation, fatigue, etc. Similar complications occur in other systems mentioned above. As a result, valid and reliable stochastic models for these systems are hard to build and maintain. They require continuous learning and adaptation based on incoming data that reflects system evolution.

The main motivation for modeling and simulating these systems is to construct good decisionmaking policies for their management. In a typical call center with multiple call types and multiple agent types (who can handle subsets of call types), one must decide how many agents to hire and train, for what call types, construct work schedules for these agents that respect union agreements, specify dynamic routing rules for arriving calls and for agents that become available, while meeting certain (stochastic) constraints on the quality of service of the systems (e.g., on the distributions of waiting times and call abandonments), and dothis at the least possible cost. Solving such stochastic optimization problems via simulation, both for long-and medium-term planning (days or months in advance) and for short-term decision making and recourse to face unexpected situations, are challenging tasks. We will illustrate these types of modeling and optimization problems with concrete examples and data, and will review some recent models and ideas.

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014

Thomas Bonald Telecom ParisTech et LINCS Réseaux de Whittle et applications Résumé : Les réseaux de Whittle sont un outil puissant et pourtant peu connu de la théorie des files d'attente, généralisant les réseaux de Jackson et les réseaux de Kelly. Ils permettent de représenter certaines dépendances entre les taux de service des différentes files d'attente du réseau tout en gardant une forme explicite de la distribution stationnaire. Cet exposé constituera une introduction à cette classe de réseaux de files d'attente. Nous en donnerons la définition et les principales propriétés, puis montrerons en quoi ils permettent d'obtenir très simplement certains résultats classiques de la théorie des files d'attente, comme l'insensibilité de la discipline de service "processeur partagé". Nous terminerons l'exposé par les applications des réseaux de Whittle aux problèmes d'ingénierie des réseaux de communication et des centres de données.

Anne Bouillard École normale supérieure Algorithms for and from network calculus Résumé : Network calculus (NC) is a theory based on the (min,plus) algebra and whose aim is to compute worst-case performance bounds in communication networks. This theory abstracts flows circulating in a network and the service offered by the networks elements by functions on which computations are performed. After a review of the basics of the theory, two problems will be addressed. First the problem of computing exact worst-case performance bounds in a feed-forward network, and second, the problem of supervision of a flow that uses some basic concept of NC but with a different aim.

Introduction

I will present a result obtained with Bruno Gaujal and Florence Perronin on simulation of grid Markov chain. Random walk are frequently used to model behavior of systems or complex randomized algorithm. one often needs to be able to compute or sample the stationary distribution. The problem we consider here is sampling the stationary distribution of a random walk on a multidimensional grid of great dimension, when its size prevents from directly computing it. A known solution is Monte-Carlo method, but it only give a asymptotically perfect approximation. Instead, I will present a method giving a perfect sample in finite time.

Model

We work on a finite grid S := {1, . . . , N} d , where both the span N and the dimension d are large, to which we add forbidden moves, couple of (point,direction) retrieved from the possibles moves on the grid. We make the hypothesis that the grid stay strongly connected. Hence there exists an unique stationary distribution.

Perfect sampling

This section is devoted to the construction of a perfect sampling algorithm of a random walk X(n) over S where certain arcs are forbidden. This is done in several steps

Coupling and Rejection

The random walk over a grid with forbidden arcs is an irreducible, finite, discrete time Markov chain over a finite state space S denoted X(n) n∈◆ , with transition matrix P. By definition, for any position a and any direction m = ±e i , P a,a+m = 1 qa where q a is the number of possible moves from a.

From (X(n)) n∈◆ , one can construct a continuous time Markov chain Y(t) t∈| over the same state space. The generator Q of Y is obtained by multiplying each line a in P by q a and defining the diagonal element Q a,a as Q aa =b q a P ab . Therefore, the rates from a to a + m are all equal to one : Q a,a+m = 1. From Y(t), it is possible to extract a new discrete time Markov chain, Y(n) n∈◆ by uniformization. Its transition matrix is Id + Λ -1 Q, where Λ (uniformization constant) is any positive real number larger than all q a 's. Since the total rate out of any state in Y is bounded by 2d, it can be uniformized by Λ = 2d.

While it can be difficult to construct a grand coupling for chain X, such a construction is easy and natural for the chain Y since the rates are all equal. To couple the walks starting from all states, just pick one direction uniformly among the 2d possibilities and make every walk move in that direction. Those for which the move is not possible stay at the same position. The stationary distribution of the first chain can be obtained from the new one through rejection. We just have to simulate one more step and reject if the new step is the same Theorem 1. If we can sample Y under its stationary distribution, then using the rejection we obtain samples distributed according to the stationary distribution of X(n).

Perfect Sampling, coupling from the past

In order to obtain a perfect sample, we need to have a sufficiently long random chain such that the walks starting from every state using this chain end at the same point, which will be the answer of the algorithm.

For that, we need to provide an efficient criterion choosing the move function of the current state and a random word. This was difficult with the first chain, but immediate with the second. If we go forward, adding moves at the end of the simulation until all trajectories converge, we add a bias caused by the dependence of coupling time on the step. Instead we need to go backward, adding the new moves at the beginning of the simulation. For every move added, the set of new trajectories is a subset of the previous one, meaning that when the last set is a point, it will stay unchanged by any addition of moves.

Lemma 1. The coupling from the past converges to a point sampled according to the stationary distribution

intervals

The size of the grid forbid us to compute every trajectories, but we can compute any supersets of the possi-ble positions. This can increase the coupling time, and could make it infinite in the case of a bad choice of superset (for example, the classical interval). In case of convergence, we know that the point obtained is the same as the one we would obtain with computing trajectories separately. We use a double superset, made of one interval, behaving as if there were no forbidden move on the grid, and a set of point, created by the forbidden moves. We can show that the expectation of the coupling time is polynomial (O(N 2 d log d)) and thus the number of point stays computable

Complexity

The goal of this section is to bound the expected time and space cost for random grid, ie a grid with a known number k of random forbidden moves. We suppose the dimension d to be smaller than the span N of the grid.

Rejection

The reject probability is the probability that the next move cannot be taken by the walker who is in a stationary state (of Y). In general this probability is bounded by 1 2d (for example if the walker has all moves forbidden but one). So that the expected number of rejections is always bounded by 2d. However this bound is very loose. In the case without any forbidden move, the expected number of tries needed is N N-1 . The numerical experiments on random graph with a number of forbidden moves smaller or comparable to N show a probability of this order.

Interval Coupling time

We have some results on random or without forbidden moves grid, as approximations. With this, we can conclude that most of the time is taken by the convergence of the interval, taking an average of O(d log dN 2 ) steps, and the average memory cost of the separated point is of the order of kN π 2 . A step can be computed in time O(k), so we obtain an algorithm for perfect sampling of grids whose time complexity is O(kd log dN 2 ).

Lemma 2 (Coalescence in dimension 1). For any

T > 0, P(C > T ) ≤ cos T π N+1 1 + O( 1 N 2 ) . Lemma 3 (Coalescence in dimension d).

Simulation parfaite dans un réseau fermé de files d'attente

Anne Bouillard, Ana Bušić, Christelle Rovetta * Inria -Département d'Informatique de l'ENS 23 avenue d'Italie, 75013 Paris, France christelle.rovetta@inria.fr

Modèle

Nous présentons une méthode par chaîne bornante pour la simulation parfaite efficace de réseaux fermés de files d'attente. On considère un réseau fermé de K files d'attente ./M/1/C à M clients. Chaque file d'attente k ∈ {1, . . . , K} := F a une capacité C k et un taux de service ν k . On note p i,j la probabilité pour un client venant d'être servi dans la file i d'être dirigé dans la file j. Une transition peut s'effectuer si la file i n'est pas vide et la file j n'est pas pleine. Le réseau est modélisé par un graphe orienté G = (F, R) avec R = {(i, j) | p i,j > 0}. On le supposera fortement connexe. On note S l'espace des états du système. Ainsi S est l'ensemble des

x = (x 1 , x 2 , . . . , x K ) ∈ N K véri- fiant K k=1 x k = M et ∀k 0 ≤ x k ≤ C k . Une borne supérieure pour |S| est donnée par K+M-1 K-1
. Soit (i, j) ∈ R, on définit t i,j : S → S la fonction qui décrit le service d'un client de la file i vers la file j :

t i,j (x) =
xe i + e j si x i > 0 et x j < C j , x sinon.

L'évolution de ce système peut être décrite par une chaîne de Markov et être échantillonnée grâce à l'algorithme 1 de simulation parfaite (voir [1]). La taille exponentielle de l'espace des états et la nonmonotonie de la chaîne rendent l'algorithme 1 inutilisable en pratique. On propose dans ce résumé une chaîne bornante sous forme de diagramme afin de pouvoir appliquer la simulation parfaite pour K et M grands. On exploite la contrainte forte du nombre de clients fixé à M, les états sont représentés dans le diagramme par des chemins partant de (0, 0) et allant à (K, M).

Représentation par un diagramme

Diagramme

On appelle diagramme complet un graphe orienté D = (N, A). Les noeuds sont placés dans une grille à (K+1) * . Auteure correspondante.

Algorithm 1: Simulation parfaite avec S Data:

(U -n = (i -n , j -n )) n∈N suite i.i.d. Result: s ∈ S suivant la distribution stationnaire begin n ← 1; t ← t U -1 ; while |t(S))| = 1 do n ← 2n; t ← t U-n • • • • • t U -1 ; end return s 0 = {t(S)} end colonnes et (M + 1) lignes. • N = {(c, ℓ), | c ∈ {1, . . . , K -1}, | ℓ ∈ {0, . . . , M}}∪ {(0, 0)} ∪ {(K, M)}. • A = {((c -1, ℓ), (c, ℓ ′ )) | 0 ≤ ℓ ′ -ℓ ≤ C c }.
On note Π(D) l'ensemble de tous les chemins de (0, 0) à (K, M). Ce sont des chemins de longueur K. On définit la fonction f : S → Π(D) qui à chaque élément x = (x 1 , . . . , x K ) ∈ S associe un chemin f(x) = ((0, 0), (1, x 1 ), . . . , (c,

φ(S) = (N, A f(S) ) et ψ(D) = f -1 (Π(D)).
Si S = ψ(D), on dit que D est un représentant de S. Si S ⊆ ψ(D), on dit que D est un sur-représentant de S.

Le diagramme complet D est un répresentant de S.

Transitions

Soit (i, j) ∈ R. On définit la fonction T i,j par :

T i,j (D) = φ • t i,j • ψ(D) pour D ⊂ D.
Proposition 1 (i) Si D est un sur-représentant de S alors T i,j (D) est un sur-représentant de t i,j (S).

(ii) Si |ψ(D)| = 1 alors |ψ • T i,j (D)| = 1.
Pour effectuer une transition T i,j sur un diagramme, on sépare ses arêtes en trois sous-ensembles :

-

Vide = {a ∈ A | a ∈ w ∈ Π(D) et f -1 (w) i = 0}, -Plein = {a ∈ A | a ∈ w ∈ Π(D) et f -1 (w) j = C j }, -Actif = {a ∈ A | a ∈ w ∈ Π(D) , f -1 (w) i > 0 et f(w) j < C j }.
Les ensembles Vide et Plein contiennent les arêtes ne pouvant pas subir la transition T i,j . L' ensemble Actif a contrario contient les arêtes pouvant subir la transition T i,j . On note Actif ′ l'ensemble des arêtes du diagramme T i,j (D) ayant subi la transition T i,j . Il est obtenu à partir de Actif en abaissant ou en remontant les arêtes des colonnes comprises entre i -1 et j -1. On obtient ensuite :

T i,j (D) = (N, A ′ ) avec A ′ = Vide ∪ Plein ∪ Actif ′ .
Une transition T i,j (D) s'effectue avec une complexité en temps en O(KM 2 ).

Simulation parfaite & diagramme

Théorème 1 Il existe une suite finie de transitions

T = T iρ,jρ • • • • • T i 1 ,j 1 telles que |Π(T (D))| = 1.
Le théorème 1 prouve que l'algorithme 2 se termine. La proposition 1 assure qu'avec une même suite (U -n ) si |Π(T (D))| = 1 alors w 0 = {Π(T (D))} = {t(S)} ce qui implique que w 0 est distribué selon la loi stationnaire. On appelle temps de couplage la valeur de n permettant de terminer l'algorithme de simulation parfaite. La figure 3 donne le temps de couplage moyen pour 200 simulations des deux algorithmes pour un même réseau. L'expérimentation avec l'algorithme 1 ne peut pas en temps raisonnable (moins de douze heures) dépasser M = 23.

U -n = (i -n , j -n )) n∈N suite i.i.d. Result: w ∈ Π(D) begin n ← 1; T ← T U -1 ; while |Π(T (D))| = 1 do n ← 2n; T ← T U-n • • • • • T U -1 ; end return w 0 = {Π(T (D))} end FIGURE 3 -K = 8 et C = M 2

Conclusion

La principale contribution est la description d'une nouvelle technique de simulation parfaite pour les réseaux fermés de files d'attente. Notre approche permet d'échantillonner la distribution stationnaire pour des réseaux ayant des files à capacité finie ou infinie, à serveurs uniques ou multiples. La démonstration de bonnes propriétés du diagramme impliquera que les transitions T i,j (D) peuvent s'effectuer en O(KM). D'autre part, les simulations montrent que les temps de couplage sont assez proches. Une avancée majeure serait de trouver une borne théorique donnant le temps de couplage de l'algorithme 2 en fonction de celui de l'algorithme 1.

Bibliographie

Introduction

Dans le milieu pétrolier, obtenir des indicateurs de fiabilité précis sur des barrières instrumentées de sécurité (HIPS) est un enjeu important. Étant donné la constante évolution du contexte opératoire de ces équipements (installation sous-marine, changement climatique, politiques de maintenances sophistiquées, . . . ), l'analyse fiabiliste de ces barrières devient de plus en plus difficile à mettre en oeuvre. De fait, les langages de modélisation usuels et les méthodes analytiques et numériques de calculs standards apparaissent de moins en moins adaptés pour prendre en compte la complexité des modèles mathématiques sous-jacents. Une solution alternative efficace est la simulation de Monte-Carlo (MC) combinée aux Réseaux de Petri (RdP) [1]. Néanmoins, obtenir des indicateurs de fiabilité précis sur ces équipements très fiables demeure un sérieux défi. En effet, la simulation d'événements rares nécessite un temps de simulation très long. Afin de répondre à cette problèmatique, des techniques d'accélération de Monte-Carlo ont été développées telles que l'Importance Sampling et le Multilevel Splitting [2] . Toutefois, ces approches nécessitent d'avoir une connaissance approfondie du modèle mathématique du système et sont fortement dépendantes du choix de paramètres qui sont difficiles à déterminer de manière automatique. Il s'ensuit que ces méthodes sont en pratique difficilement adaptables aux RdP. Une technique moins courante consiste à tronquer à droite les distributions des événements menant le plus directement à l'événement d'intérêt. On se réfère à la Méthode de Conditionnement Temporel [3] (MCT), également appelée Failure Forcing [4] dans le cadre markovien.

Cette technique ne nécessite pas de connaître les lois intervenant dans le modèle, mais est uniquement définie lorsque l'événement d'intérêt est absorbant. Par conséquent, afin de remédier à ce problème, nous introduisons dans un premier temps une extension de MCT (EMCT) pour des modèles simples représentant des cycles de vie répétés où la défaillance qui mène à l'événement rare est soit directe, soit en conflit avec d'autres types de pannes. Ensuite, nous proposons une procédure informatique, valable uniquement si les composants sont indépendants, et nous la combinons enfin à la EMCT. Par suite, notre but est d'évaluer la Probability of Failure on Demand (PFD) de ces HIPS quelle que soit leur structure. Ainsi, nous présentons des résultats numériques qui illustrent le potentiel de ces méthodes.

Analyse fiabiliste

Dans le domaine de la sûreté de fonctionnement, l'une des principales quantités d'intérêt mesurant les performances d'un HIPS est la PFD. Nous proposons la définition suivante. Définition 2.1 L'indicateur PFD(t) correspond à l'indisponibilité instantanée à l'instant t d'un système relatif à la sécurité le rendant incapable d'accomplir correctement sa fonction de sécurité. Cette quantité est définie de la manière suivante pour tout t ∈ R + :

PFD(t) = P (Φ(X(t)) = 0) où :
-

X(t) = (X 1 (t)
L( ξ(k) ) = L ξ (k) | min ξ (1) , . . . , ξ (p) ≤ ∆ .
Par suite, en appliquant le Principe 3.1, nous avons calculé les PFD associées aux systèmes de la Section 2, ce qui nous a permis d'obtenir un facteur de déconditionnement pour chaque structure.

Méthode de Dissociation

Cette procédure qui peut être associée à la méthode de Monte-Carlo (MC) classique ou combinée à EMCT est décrite de la manière suivante. Principe 4.1 Nous supposons les m composants du système indépendants. Au lieu d'observer le comportement conjoint de ces m composants au cours de n simulations, cette méthode, dite de Dissociation, consiste à observer chacun des m comportements individuellement. Ainsi, la première étape revient à simuler n fois chaque processus univarié indépendamment. Dans un deuxième temps, nous combinons chacun des comportements individuels entre eux afin d'obtenir le comportement global du système. En d'autres termes, la seconde étape est une phase de posttraitement des n simulations, qui permet de croiser tous les cas possibles afin d'obtenir au final n m observations. Cette méthode présente un intérêt si nous considérons, par exemple, un système de deux composants réparables avec deux réparateurs disponibles. Cependant, dans le cas où un seul réparateur est à disposition, elle ne pourra pas être appliquée.

Exemples numériques

Pour chacune des structures présentées dans la Section 2, nous avons estimé la PFD à l'instant t = 8759 de systèmes constitués de chaque composant, la loi de défaillance de la 1 ère panne suit une loi Exp(5 × 10 -6 )2 et de la 2 ème une loi Exp(1.74 × 10 -4 ). Les deux pannes sont réparées selon une loi Exp(1.742 × 10 -1 ). Pour chacune des méthodes, l'estimateur est obtenu après 1 000 000 de simulations. Lorsque l'événement n'est pas très rare (PFD(t) > 10 -5 ), réaliser 1 000 000 de simulations semble suffisant pour observer que EMCT donne de meilleurs résultats que MC pour un temps de calcul équivalent. Par contre, lorsque la quantité d'intérêt devient très rare (PFD(t) ≤ 10 -9 ), ce nombre de simulations ne permet pas d'obtenir d'estimateurs avec MC et EMCT. Cependant, nous pouvons observer que les combinaisons avec la méthode de Dissociation permettent d'améliorer significativement les résultats obtenus quel que soit le degré de rareté de l'événement d'intérêt. Des algorithmes sont en cours de traitement afin de moyenner diffèrentes estimations pour obtenir des estimateurs plus stables.

Introduction

L'évaluation des performances de très grands réseaux et systèmes nécessite parfois une simulation très longue de modèles markoviens, que ce soit par exemple pour capturer des séquences d'événements rares, ou bien pour évaluer des taux de couverture, ou encore pour analyser des basculements de comportement dans des modèles "raides". Pour réduire la durée de la simulation d'une seule trajectoire longue, l'idée est d'utiliser plus efficacement les ressources à disposition, en particulier les ressources de calcul (coeurs), en anticipant les calculs sur des tranches de temps se recouvrant au minimum. L'argument clé est, sous certaines conditions, l'efficacité du couplage de trajectoires qui permet de garantir l'obtention de la trajectoire complète. L'objectif de ce travail est donc de proposer un algorithme parallèle de simulation de trajectoires longues basé sur une technique de "couplage horizontal" de trajectoires, d'analyser son coût et de donner les premiers résultats d'efficacité.

Simulation de modèles markoviens

On considère un modèle markovien à temps discret {X t } t∈◆ décrit par sa fonction de transition Φ : X t = Φ(X t-1 , e t ) où e t est l'événement se produisant à l'instant t (généré aléatoirement). La simulation classique d'une trajectoire de longueur T consiste, pour chaque instant, à générer aléatoirement l'événement selon les probabilités du modèle puis à mettre à jour l'état selon la fonction de transition (algorithme 1). Une technique de parallélisation consiste à calculer séparément des segments disjoints de la trajectoire, en anticipant l'état de départ dans chacun de ces segments. Cette technique initialement proposée par Nicol et al. [4] ou Fujimoto [3] est la suivante : La période simulée est partitionnée en P tranches de * . Ce travail a été partiellement financé par l'ANR Marmote Algorithme 1 : Algorithme séquentiel Données : x 0 (état initial), T pour t de 1 à T faire e t ← genere_evenement(); x t ← Φ(x t-1 , e t ); temps de durée fixe ∆ = T/P démarrant à t i = i∆ pour 0 i P -1 affectées respectivement aux P coeurs. Pour chaque coeur on génère également la séquence des événements et l'état initial à l'instant t i . Plusieurs phases s'ensuivent durant lesquelles : chaque coeur i simule la trajectoire du segment [t i , t i+1 ] selon l'algorithme 1 en partant de l'état x t i ; les états de début de tranche sont corrigés par la valeur obtenue en fin de la tranche précédente à la fin de la phase précédente. Les phases sont répétées jusqu'à ce que les états de début de tranche ne soient plus modifiés. Il est clair que dans le pire des cas cet algorithme simulation parallèle nécessite autant d'étapes de calcul que l'algorithme séquentiel (au pire P phases). Le temps de calcul peut encore être réduit au prix de la précision en simulant un ensemble d'états possibles au lieu d'une trajectoire exacte [2]. L'idée est de remplacer le calcul de l'état par le calcul de bornes sur cet état : ainsi, les conditions aux limites des tranches de temps sont moins strictes (inclusion) ce qui réduit le nombre de phases de l'algorithme. Cela permet d'obtenir rapidement une estimation de la partie de l'espace d'état atteinte par le modèle.

L'algorithme Catch me if you can

L'algorithme que nous proposons vise à adapter la méthode développée par Nicol et al. en conditionnant l'arrêt du calcul de la trajectoire sur une tranche de temps par une condition de couplage.

Couplage de chaînes de Markov

Considérons une séquence d'événements {e 1 , e 2 , • • • , e T } qui permet la génération de deux trajectoires {X a t } et {X b t } par applications successives de Φ à partir des états initiaux x a 0 et x b 0 .
Étant donné la nature du schéma itératif, si les deux trajectoires couplent à l'instant t alors elles seront confondues après t. Le premier instant de rencontre des deux trajectoires est appelé instant de couplage, noté τ.

Principe de parallélisation

On suppose que la séquence des événements est fixée et partagée en lecture par tous les coeurs. Elle peut être calculée au préalable ou à la volée en parallèle. À chaque coeur i, 0 i P -1, on attribue un intervalle de temps de taille variable pour calculer un fragment de trajectoire. C'est à dire une date de départ et un état de départ x t i correspondant (fixé arbitrairement ou par une heuristique). Chaque coeur i calcule donc en parallèle la trajectoire issue de x t i à partir de la date t i . Le calcul est arrêté dès que la trajectoire couple avec une trajectoire déjà calculée par un autre coeur, ce qui permet d'éviter des calculs redondants. C'est pourquoi cet algorithme s'appelle catch me if you can. L'algorithme proposé est illustré en Figure 1. Si le principe de cet algorithme est simple, les difficultés résident dans sa mise en oeuvre. En effet, même si l'on travaille en mémoire partagée, le test compare des données associées à des coeurs différents, ce qui correspond à des synchronisations en mémoire qui peuvent dégrader fortement les performances de l'algorithme. Du point de vue théorique, si on note N i le nombre d'itérations effectuées par le coeur i avant couplage, la durée totale de simulation (makespan) en nombre d'itérations est N = max i N i . Pour des t i+1t i de l'ordre de T P on peut montrer que :

N ≃ T P + α,
où α correspond au surcoût lié au couplage de trajectoires et peut être majoré stochastiquement par le maximum de P temps de couplage, ce qui croît très lentement en fonction de P et laisse espérer au moins une efficacité asymptotique.

Premières analyses

L'un des critères de performance de l'algorithme est le nombre d'itérations (calculs de la fonction Φ) nécessaires à l'obtention de la trajectoire complète. L'exemple de la figure 2 illustre l'accélération obtenue en parallélisant le calcul de la trajectoire. Dans ce modèle le temps moyen de couplage est grossièrement borné par une fonction linéaire en la capacité de la file [1] (ici de l'ordre de 250 itérations), ce qui est négligeable devant T P . Par contre si l'on augmente le nombre de coeurs le temps de couplage devient non-négligeable. Cet effet sera accentué par le coût lié aux contraintes de synchronisations entre les différents coeurs. Dans notre cas de simulation un nombre de coeurs de l'ordre de 10 donnerait un bon compromis entre accélération et utilisation de ressources.

Analyse de différents types de systèmes

Perspectives

La technique d'agrégation temporelle que nous avons implémentée dans Ocelotl permet de représenter une grande quantité d'information sous la forme d'une description macroscopique simple dans des délais adaptés à une analyse interactive. Nous aimerions, cependant, pouvoir représenter conjointement les ressources, en particulier dans les cas hétérogènes, afin de déterminer leur influence sur les performances de l'application. Nous sommes ainsi focalisés actuellement sur une extension de ces travaux à l'agrégation simultanée de l'espace des ressources et du temps.

Bibliographie

Introduction

With the diversity of applications and technologies present in communication networks, the study of their performance has become a complex issue. The large number of parameters to take into account to model communication network performance increases the risk of spurious associations between the explanatory variables. Therefore studies based on correlation become a challenging approach. On the other hand, causal models, and their representation as Directed Acyclic Graphs (DAGs), offer simple and attractive models. Using simple graphical criteria, one can predict interventions on the system by using data collected through passive measurements. In this paper, we present an example of a causal study of communication network performance. We place ourselves in an emulated environment, where our predictions can be verified, and we show the benefits of the approach.

Causal study of an emulated network

For this paper we study the performance of TCP observing FTP traffic. We focus on the network performance, represented by the TCP throughput. To be able to verify our predictions we emulate a network using the Mininet software [1]. The experiment consists in one client that downloads files from a FTP server, both machines are connected to different routers (R1 and R2 respectively). In order to create a more realistic scenario we add two other machines, one connected to R1 and the other to R2, creating cross traffic. The FTP traffic is recorded at the server side using Tcpdump tool and the dataset that we obtain is presented in Table 1.

Using the PC algorithm [4] with the Hilbert Schmidt Independence Criterion (HSIC) [5], we obtain, from the data summarized in Table 1, the causal model represented Figure 1. While TCP is a feedback based protocol, it is important to notice that we are studying the causal relationships between the parameters defining our system. In this case we obtain a causal graph, rep-FIGURE 1 -Causal model of emulated FTP traffic resented as a DAG, where the causes of the throughput are disclosed and its direct causes are represented as its parents. In this model we observe the receiver window (rwin), delay (rtt), and loss (retrscore), the empirical causes of the throughput, that are, indeed, found as its direct parents. Due to space constraints we cannot describe all the properties of this model and focus on prediction. We want to estimate the distribution of the throughput after an intervention on the retrscore, where its value is set to 1%.

To estimate the total effect, on a parameter Y, of an intervention on a parameter X, we use the Back-door Criterion [3] defined as follows.

Definition 1 (Back-door criterion) A (set of) variable(s) Z is said to satisfy the back-door criterion relative to an ordered pair of variables (X i ,X j ) in a DAG G if : (i) no node in Z is a descendant of X i ; and (ii) Z blocks every path between X i and X j that contains an arrow into X i .

Using this definition, we can predict the effect of an intervention as follows.

Theorem 1 (Back-door adjustment) If a set of variables Z satisfies the back-door criterion relative to (X, Y), then the causal effect of X on Y is identifiable and given by the formula :

P(y | do(X = x)) = z P(Y = y | X = x, Z = z)P(Z = z).
where the do() operator represents the intervention of setting the parameter X to the value x We estimate the probability density functions (pdfs) by first estimating the marginals using normal kernels and then modeling the multivariate and conditional pdfs using T-Copulae [2].

To predict the effect of an intervention on the retrscore parameter, using the causal model presented in Figure 1, we apply the Back-door adjustment formula (Theorem 1) with the blocking set Z = {nlac, rtt}. To verify our prediction we set up a new experiment where we drop 1% of the packets at the router R1 and compare the throughput we obtain with the one that was predicted using the equation from Theorem 1. Figure 2 presents the probability density functions corresponding to the throughput prior to intervention (in dash-dot line), the post-intervention throughput estimated with the Back-door criterion (dashed line) and the throughput observed when we manually modify the loss (solid line). While not perfectly matching, the graph-based prediction (dashed line) and the real density after manual intervention (solid line) show similar shapes. In particular, both densities have a single mode around the same value (∼ 100 kbps). The fact that the experimental throughput values after intervention (solid line) show less variance can be explained by the small number of samples (20) used to make the estimation. Also our method to estimate the post-intervention throughput (dashed line) uses normal kernels and a T-copula which tend to widen and smooth the estimated post-intervention density.

Concluding remarks

Due to space constraints, many aspects of this study could not be presented. The study of the causal model we obtain discloses many properties of the system that would need further explanations and would lead to future studies. However it should be noticed, from the Back-door adjustment equation (Theorem 1) that the values present in our original dataset dictate the range of predictions that can be made. We are now working on adopting parametric modeling of the probabilities of the parameters that would allow to overcome this limitation. Most traffic models used to estimate the performance of a cache server are based on a fixed document catalog. As an example, the Independent Reference Model (IRM) further assumes that all requests are i.i.d., which allows to calculate theoretical estimates of the hit ratio, defined as the number of cache hits over number of requests. However, in many applications, content is dynamic: It appears at a certain instant, its popularity varies over time and can eventually disappear. In this paper, we aim at building a tractable model that captures this feature (see [3] for details).

Semi-experimental driven modeling

We base our analysis on two datasets. The first dataset, called YT, comes from three months of YouTube traffic of the Orange Network in Tunisia. The second, called VoD, collects over three years requests from the Orange VoD service in France. To discover the relevant features of our request sequences, we apply the semi-experimental methodology [2]. These semi-experiments aim to break specific structures of the request process by means of random shuffling; we then compare the resulting hit ratio, computed via simulation, with that of the original sequence and determine the importance of the broken structure. They consist in: 1) shuffling all requests, thus breaking all correlations; 2) shifting each document request sequence randomly, thus breaking the correlation between first requests; 3) fixing the first and last request of each document and shuffle all requests in between, thus breaking correlations within this sequence. The results of semi-experiment (1) showed that the hit ratio curve differs considerably from that of the original sequence, thus confirming that IRM does not capture all the features of our datasets. In the case of ( 2) and ( 3 where the random variable τ d denotes the document lifespan. We further assume that the distribution of the pair (λ d , τ d ) is stationary. We consolidate the two layers into the total request process, which is the superposition of the document request processes. Figure 2 shows a schematic view of our model.

Hit ratio estimates and model validation

In order to estimate the hit ratio, we introduce the auxiliary process X t , that counts the number of different requested objects up to time t. We then observe that for a LRU cache of size C, if an object is requested at time zero and later at time Θ > 0, the latter request will be a hit if and only if X Θ < C. Equivalently, the request at time Θ will be a hit if and only if Θ < T C , where T C is the first passage time of process X at level C. We then invoke the so-called "Che approximation" [1], in which we assume that the random variable T C is concentrated and thus can be approximated by a constant t C . To calculate the constant t C , we notice that, by definition of T C , we have the equality X T C = C. Thus, we impose that the constant t C must satisfy this equality in mean, that is, E [X t C ] = C. Process X is a non-homogeneous Poisson process. Define its mean function by Ξ(t) = E [X t ]; as Ξ is an increasing function, t C is then given by the formula t C = Ξ -1 (C). We make Ξ explicit in the following proposition.

Proposition 1

The mean of the process X is given by

Ξ(t) = γ E 2t + 1 -e -λt τ -t - 2 λ 1 {τ≥t} + γ E 2τ + 1 -e -λτ t -τ - 2 λ 1 {τ<t}
where (λ, τ) is distributed as any (λ d , τ d ).

H d =      λ d τ d -1 + e -λ d τ d (λ d τ d -1)(1 -e -λ d t C ) + λ d t C e -λ d t C if τ d < t C and τ d ≥ t C , respectively.
For each document, unbiased estimates for λ and τ can be easily calculated; call these estimators λ and τ, respectively. Such estimators are valid only for documents with a number of requests n ≥ 2; we cannot neglect, however, the documents that have only one request, since they compose a considerable part of the data. Due the fact that they add only misses, we can still incorporate them into the hit ratio as follows:

HR = E [H d ] E [n d ] = E [H d | n d ≥ 2] E [n d | n d ≥ 2] + P (n d = 1) P (n d ≥ 2)
.

Using the previous estimates, we compute t C and

E [H d | n d ≥ 2
] by plugging estimates τ and λ into the formulas obtained in Proposition 1 and 2 and take averages.

We observe the result of this estimation on the VoD dataset in Figure 3, where we compare it to the actual hit ratio of the request sequence, obtained via simulation, and the estimation obtained via the "classic" Che approximation under the IRM setting. We observe that the hit ratio is noticeably underestimated in the latter case, whereas our model fits well the real hit ratio.

Introduction

Content distribution and in-network caching have emerged as key features of several network architectures to accommodate the current content usage patterns (Video on Demand, User-Generated Contents, Web browsing, the recent Dynamic Page Caching concept of Akamaï, etc.) while reducing congestion and improving access speed as networks increase in size. In order to guarantee the latter performance, network designers and engineers need tools to quickly gain insights on the behaviour of the multi-cache systems that arise from their contentoriented architectures when deploying caches across the networks.

The analysis of these cache networks is significantly challenging due to their irregular topologies, the heterogeneity of their nodes (different replacement policies and capacities), and the statistical correlation of requests. Few modelling attempts have been proposed to solve this problem on general and heterogeneous Time-To-Live (TTL)-based cache networks where requests are described by Markovrenewal processes [1] or Markov-arrival processes (MAPs) [2]. However, the latter derivations suffer from a large complexity which limits their practical interest.

In this work, we aim at presenting an analytical methodology to approximate the performance metrics of heterogeneous networks where (i) caches are running the Least Recently Used (LRU), First-In First-Out (FIFO), or Random (RND) replacement algorithms, (ii) all requests are routed as a feed-forward (or hierarchical) network, (iii) streams of requests are described by simple MAPs [3] as briefly explained in Section 2.1.

Single cache approximation (SCA)

We consider a cache system that can accommodate C files requested from a catalog of size K (K > C).

We assume that the cache is connected to a server where all files are permanently stored such that requests arrived first on the cache and the missed ones are forwarded to the origin server. We also assume a zero delay cache-server system i.e. request processing/forwarding and file downloading times are negligible in comparison to the inter-request times [1].

In this section we present a simplified traffic model which accounts for request correlations, then we recall the TTL-based model of LRU and FIFO/RND caches, and we provide the approximate metrics of interest (hit and occupancy ratios).

Workload model

Let us denote by {R k , k = 1, . . . , K} the general point process describing the sequence of requests for file k on the cache. Here, we aim at providing a simple and accurate approximation of R k based on a minimal available information that engineers could easily measure or estimate at edge nodes of the network.

As inputs of our process model, we choose the request rate, the squared coefficient of variation (scv), the skewness factor, and the lag-1 (and possibly lag-2) autocorrolation. These quantities are general enough to capture the main statistical parameters (per-file popularity, temporal locality, and correlation coefficients) of the arrival process R k with a MAP (Cf. the MAP-Match procedure in [3]). Therefore, all request streams in this paper will be approximated by the switched MAP(D 0 , D 1 ) described in [3,Sect.4].

Cache model

We recall a result from [1] that the characteristic time (i.e. the maximum inter-request time of a given file that yields a cache hit) converges in distribution (under some conditions that hold in practice) to deterministic and exponentially distributed random variables as K ↑ ∞ for LRU/FIFO and RND caches respectively under general stationary request processes. Hence LRU, FIFO and RND caches may be analysed through their corresponding TTL-based models. In this work, we consider that TTLs are i.i.d with an Erlang distribution E r (m, µ) (or a PHdistribution PH(α, T) with m transient states). For RND and LRU/FIFO caches we set m = 1 and m ≈ 10 ≫ 1 (sufficiently large to approximate a constant) respectively.

Performance metrics

In this section, we focus on the calculation of the file hit {H k } K k=1 (resp. occupancy {O k } K k=1 ) probabilities defined as the stationary probability that a request for file k yields a cache hit (resp. the stationary probability that file k is in the cache at any time). To do so, we should find the rate µ of the TTL distribution. This is done by solving the following equation

K k=1 O k = C.
Based on the exact results derived in [1,2] we are able to derive the closed-form expressions of H k (see [2,Lemma 5]) and O k (see [1, Prop. 2.9]) needed in the latter equation.

Networks with feed-forward routing

In this section, we generalize previous results to cache networks where requests flow in the same direction (i.e. a cache does not query a cache from which it receives missed requests) by describing the main network traffic transformations in cache network deployments where requests are correlated and routed on a unique feed-forward network.

Miss process characterization

The SCA extends easily to tandem or linear cache networks. Indeed, one need to characterize the miss process of each cache, then apply the MAPMatch procedure on the miss streams, and calculate the metrics of interest as done in Section 2. If request processes are MAPs, so are miss processes of PH-distributed TTL models [2, Thms 3 & 4]. However, the complexity of the exact miss process characterization [2] is significantly reduced with our workload model.

Multiple sources

A cache may receive requests from several (not necessarily independent) sources such as caches or users. This is the common case for tree networks. The exact characterization of the aggregation of N MAPs requires a strong independence assumption among request streams being superposed, an important calculation to evaluate the Kronecker sum, and a huge memory to store the resulting MAP (limitations of the exact method in [2]). In this work, we propose to calculate the parameters needed for MAP matching [3] as weighted sum of corresponding parameters of MAP components (e.g the resulting scv is

c 2 v = N n=1 λn Λ c 2 n
where λ n and c 2 n are the rate and the scv of the n-th MAP being aggregated, Λ = N n=1 λ n is the total request rate). Using the MAP matched of the overall process and the SCA, we obtain the performance metrics.

Several destinations

In many situations, a cache may select the next hop among two or more caches to forward its missed requests. The latter task may be performed based on some available informations (e.g. fault tolerance, multiple shortest paths, load balancing, etc.). We denote by r j the probability that the outgoing-link j is chosen for request forwarding such that j r j = 1 (e.g. r j may be the fraction of time the outgoing-link j is up while others are down). If the cache miss processes are MAPs, so are the request processes on each link j (see [2,Lemma 4]). Finally, the MAPMatch procedure [3] is applied on each r j -thinned request process.

Conclusion

In this paper, we proposed a methodology that can be used to quickly estimate the performance metrics of single TTL-based (e.g. LRU, FIFO, RND) cache when request streams are approximated by special Markov Arrival Processes. We have also explained how our modelling attempt extends heterogeneous cache networks where all requests are routed on a unique/same feed-forward network (e.g. linear, tree, polytree) built on top of the network topology. Ongoing work are devoted to provide a detailed description of the latter extensions, to perform an extensive evaluation of our models, and to investigate the case of general and heterogeneous caches networks with arbitrary routing.

Introduction

Delay tolerant networking (DTN) [1] is an architecture that promotes data transmission even in absence of end-to-end connectivity at a given time. In DTNs, forwarding of a message is performed through mobile nodes in an opportunistic manner. To cope with intermittent connectivity and to increase delivery rate, DTN routing protocols commonly appeal to multi-copy forwarding, that is the message is replicated to relay nodes and it is delivered if one of the copies reaches the destination. The two-hop routing scheme [2], which we consider in this paper, allows the source to forward copies of the message to any idle node in contact, and a relay that receives the message from the source can transmit it only to the destination. However, in DTN applications, when mobile devices are controlled by selfish owners, DTN nodes can be unwilling to spend their own resources serving as relays or can try to gain benefits from cooperation. Some form of virtual currency is proposed by creditbased incentive mechanisms [3] to encourage communication among different nodes in DTNs. The studies in this work are centred on such type of incentives, namely addressing the reward-based incentive scheme introduced in [4] for two-hop routing. Each node that accepts the message has a certain probability to reach the destination, and the incentive scheme proposes to reward the relay that is the first one to deliver. The success probability and expected reward estimated by a relay depend on certain informations that the source announces to a relay in contact. This paper investigates the impact of the announced information on the expected reward the source has to pay, and the aim is to find a reward policy for the source that is optimal in costs. The results are obtained in two cases. First, for two relays and a decreasing density function of intercontact times between relays and the source (destina-tion), it is shown that the optimal reward policy corresponds with adapting the information the source can give to the relays. This adaptive strategy is shown to be of threshold type. Second, for an arbitrary number of relays and exponentially distributed inter-contact times, it is proved that the expected reward the source has to pay remains the same within fixed setting of informing the relays. The source can, however, reduce the expected reward, by adapting the announced information according to the meeting times with the relays.

System Model and Problem Description

The network considered consists of one fixed source node, one fixed destination node and a finite number of relays. The relays move according to a mobility model with i.i.d. inter-contact times between a relay and the source (resp. destination) corresponding to density function f s (resp. f d ). The source generates a message at time 0 and attempts to forward it to the destination through the relay nodes by replicating the message, under the twohop routing scheme. Forwarding of the message imposes a certain set of costs for a relay node :

1. reception cost, C r , is a fixed energy cost for receiving the message from the source, 2. storage cost, C s , is the cost incurred per unit time a relay carries the message, 3. delivery cost, C d , is the fixed cost for transmitting the message to the destination As an incentive for cooperation, the source proposes to a mobile node a reward that promises to compensate costs that the node expects to incur if it accepts the message. The source is assumed to be not informed whether the message is delivered by any relay, and continues to offer message copies to relays when encountering. The expected costs and the expected reward estimated by a relay depend on the information the source proposes to the relays. A mobile node encountering the source knows the time when the transmission process started. In addition, the node may be given the history of forwarding process according one of three settings that the source can use :

1. full information, when each mobile node is informed about the number of already existing copies and the corresponding contact times, 2. partial information, when each mobile node knows only the number of existing copies, 3. no information, when the relays know only their own meeting times with the source.

The reward a relay gets is a random variable depending on the stochastic process describing contacts with the destination. Moreover, the information the relay is given by the source may distinctly affect its estimation of the expected reward. Thus, within each of three settings and based on the node mobility model, the source needs to determine an amount of reward it should propose to a mobile node when encountering it, and then deploy an optimal policy of informing the nodes in order to reduce the expected reward it should pay for message delivery.

Two Relays, Decreasing Function of Inter-Contact Times

Consider a case of two mobile nodes. Assume a decreasing density function of inter-contact times between the relays and the source (destination). For fixed s 1 and s 2 , denote by s the vector (s 1 , s 2 ).

Expected Cost for a Relay

Define V (k) i (s) as the net cost for relay i = 1, 2 under the setting k ∈ {F, P, N}, and let R (k) i (s) be the reward asked by this relay to the source under this setting. Then, the expected cost for relay i takes the form,

E[V (k) i (s)] = C r + C s E[ Td ] + [C d -R (k) i (s)]p (k) i (s),
with Td being the residual inter-contact time for a relay to meet the destination, and p (k) i (s) denotes success probability estimated by relay i under the setting k when it meets the source.

Expected Reward Paid by the Source

The expected reward paid by the source for given contact times s = (s 1 , s 2 ) and under setting k, is as follows

R (k) S (s) = R (k) 1 (1 -p 2 (s)) + R (k) 2 p 2 (s),
where p 2 (s) is the probability of success of the second relay given s 1 and s 2 .

Adaptive Strategy

It is always beneficial for the source to give information to the first relay independently of s 1 .

Proposition 1 R

(F)

1 (s) = R (P) 1 (s) ≤ R (N)
1 (s). The next proposition says that between the choice of informing a relay that is the second one and not giving this information, it is better for the source not to give this information.

Proposition 2 R (N) 2 (s) ≤ R (P) 2 (s).
Comparing the settings of full information with that of no information, the following proposition shows that for the source the adaptive strategy is of the threshold type. Define the difference of the success probabilities as a function of s 1 and s 2 ,

g(s 1 , s 2 ) = p (N) 2 (s 1 , s 2 ) -p (F) 2 (s 1 , s 2 ), (1)
then for the source, it will be better to give information when g(s 1 , s 2 ) < 0.

Theorem 1 There exists

0 ≤ θ 1 < ∞ such that 1. if 0 ≤ s 1 < θ 1 , then g(s 1 , s 1 + v) ≥ 0, ∀v ≥ 0 ; 2. if θ 1 < s 1 < ∞, then (a) g(s 1 , s 2 ) < 0, ∀s 2 ∈ [s 1 , s 1 + ω(s 1 )), (b) g(s 1 , s 2 ) > 0, ∀s 2 ∈ (s 1 + ω(s 1 ), ∞),
where θ 1 is a solution of the equation g(s 1 , s 1 ) = 0 and ω(s 1 ) is a solution of g(s 1 , s 1 + v) with respect to v when g(s 1 , s 1 ) < 0.

N Relays, Exponential Inter-Contact Times

Assume the inter-contact time between a relay and the source (resp. destination) is exponentially distributed with rate λ (resp. µ).

C 2 = C r + C s /µ.
The average reward that the source has to pay is, thus, the same if the source adheres to a fixed informing mode and depends only on the number of involved mobile nodes. This reward can be, however, reduced by adapting the announced information according to meeting times with the nodes.

Introduction

Abandonment or reneging takes place when customers, unsatisfied of their long waiting time, decide to voluntarily leave the system. It has a huge impact in various real life applications such as the Internet or call centers, where customers may abandon while waiting in the queue. In this paper we study the phenomena of abandonments in a single server queue with multiple classes of users with the objective of finding resource allocation strategies. In the presence of abandonments and/or convex holding cost, a characterization of the optimal control is out of reach. We therefore develop approximations to tackle the problem, in particular, we study Whittle's index. Index Rules have enjoyed a great popularity, since a complex control problem whose solution might, a priori, depend on the entire state space turns out to have a strikingly simple structure. In a seminal work, Whittle introduced the so-called Restless Multi Armed Bandit Problems (RMABP), see [2]. In a RMABP all bandits (in our study a bandit is a class of customers) in the system incur a cost, the scheduler selects one bandit to be made active, but all bandits might evolve over time according to a stochastic kernel that depends on whether the bandit was active or frozen. The objective is to determine the control policy that, based on the entire state-space description, selects the bandit with the objective of optimizing the average performance criterion. Whittle introduced an approximate control policy of index-type, which is nowadays referred as Whittle's index. The multi-class single server queue with abandonments can be modeled as a RMABP problem and Whittle's index can be derived. The reminder of the paper is structured as follows : in Section 2 the model under study is described, and the objective of the study is presented, in Section 3 the Whittle index is derived and in Section 4 the performance of this index policy is measured. For additional results and proofs we refer to [1].

Model Description

We consider a multi-class single-server queue with K classes of customers. Class-k customers arrive according to a Poisson process with rate λ k and have an exponentially distributed service requirement with mean 1/µ k , k = 1, . . . , K. We denote by ρ k := λ k /µ k the traffic load of class k, and by ρ := K k=1 ρ k the total load to the system. We model abandonments of customers in the following way : any class-k customer in the system that has not completed service, abandons after an exponentially distributed amount of time with mean 1/θ k , k = 1, . . . , K, with θ k > 0. The server has capacity 1 and can serve at most one customer at a time, where the service can be preemptive. At each moment in time, a policy ϕ decides which class is served. Because of the Markov property, we can focus on policies that only base their decisions on the current number of customers present in the various classes. For a given policy ϕ, N ϕ k (t) denotes the number of class-k customers in the system at time t, (hence, including the one in service), and

N ϕ (t) = (N ϕ 1 (t), . . . , N ϕ K (t)). Let S ϕ k ( N ϕ (t)) ∈ {0, 1} represent
the service capacity devoted to class-k customers at time t under policy ϕ. The constraint on the service amount devoted to each class is

S ϕ k ( n) = 0 if n k = 0 and K k=1 S ϕ k ( n) ≤ 1.
The above describes a birth and death process. Let C k (n k , a) denote the cost per unit of time when there are n k class-k customers in the system when class k is not served (if a = 0), or when class k is served (if a = 1). We further introduce a cost d k for every class-k customer that abandoned the system when not being served. We then want to find the optimal scheduling policy ϕ under the average-cost criteria, that is, lim sup

T →∞ K k=1 1 T E T 0 Ck ( N ϕ (t), S ϕ k ( N ϕ (t))) dt , (1) where
Ck (n k , a) := C k (n k , a) + d k θ k n k in state n k
and under action a.

Whittle's index

The approach by Whittle is based on relaxing the original problem, allowing one bandit being made active on average, that is, lim sup

T →∞ 1 T T 0 K k=1 S ϕ k ( N ϕ (t))dt ≤ 1. (2)
This allows to decompose the original control problem into individual problems for each bandit. Whittle's index can then be interpreted as the Lagrange multiplier of the constraint such that a given state joins the passive set. The objective is now to determine the policy that solves (1) under Constraint (2). This can be solved by considering the unconstrained control problem lim sup

T →∞ K k=1 1 T E T 0 Ck ( N ϕ (t), S ϕ k ( N ϕ (t))) -W(1 -S ϕ k ( N ϕ (t))) dt , (3)
where W is the Lagrange multiplier and C(•, •) is assumed to be convex and increasing. We observe that the multiplier W can be interpreted as a subsidy for passivity.

In summary, the relaxed optimization problem can be written as K independent one-dimensional Markov decision problems, namely :

lim sup T →∞ 1 T E T 0 C(N ϕ (t), S ϕ (N ϕ (t))) -W(1 -S ϕ (N ϕ (t))) dt . (4)
Under an indexability property (established in [1] for the present model), one can define class-based indices (the so-called Whittle's index) such that the solution to ( 1) is obtained by activating in every moment those classes whose current index is larger than the Lagrange multiplier. We further prove that the optimal policy is monotone, it is fully characterized by a threshold n such that the passive action is prescribed for states m ≤ n, and the active for states m > n. We will then denote this policy as ϕ = n. We can now state the main result.

Proposition 1

The Whittle index for problem (4) is

W(n) = E( C(N n , S n (N n ))) -E( C(N n-1 , S n-1 (N n-1 ))) n m=0 π n (m) -n-1 m=0 π n-1 (m) (5)
with π n (m) the steady-state probability under threshold policy n.

The heuristic for the original K-dimensional problem known as Whittle's index policy, prescribes to serve the class k having currently the highest non-negative Whittle's index W k (n k ).

We now present Whittle's index for linear cost.

Proposition 2 Assume linear holding cost C k (n k , a) = c k n k . Then, Whittle's index for class k, with ck = c k + d k θ k , is W k (n k ) = ck µ k θ k , for all n k . (6)
W(n) W ∞ (n) w(n) Gcµ (21) FIGURE 1 -Performance of Whittle's index policy.
Whittle's index for a more general setting (with convex holding cost) as well as limiting properties are analyzed in [1]. This study allows to recover indices that have been obtained in the literature. Observe that ( 6) is independent on the number of class-k customers, which is not the case for convex holding costs.

Numerical Results

In this section we numerically analyze the performance of Whittle's index policy. In Figure 1 we plot the relative sub optimality gap of the W(n) index as well as other indices obtained in [1]. The optimal policy is computed using a Value Iteration algorithm. We carried out extensive simulations and were able to conclude that Whittle's Index policy is optimal as the workload increases, whereas indices that have been computed for multi-class systems without abandonments behave poorly.

Introduction

In the last decades, computers, operating systems and software running on it have reached such a level of complexity that it has become very difficult (not to say impossible) to control them perfectly and to know every detail about their configuration and operation mode. As a consequence, it becomes less and less reasonable to consider computer systems as deterministic, especially when measuring execution times or performance of large distributed computing systems. These systems have become so complex that it is almost impossible to fully understand their behavior, making the full reproduction of measurements extremely difficult. To some extent, although these systems have been designed and built by humans, studying modern computers has become very similar to studying a natural phenomenon. Just like most other scientific domain, many of the conclusions are based on experiment results and their analysis. Nevertheless, in many computer science articles, the experimental protocol is rarely detailed enough to allow others to reproduce the study and possibly build upon it. This is all the more surprising as reproducibility and falsifiability are yet the basis of modern science as defined by Popper. In the recent years, the more and more frequent discovery of frauds or mistakes in published results has shed the light on the importance of reproducible research in computational sciences. Many tools partially addressing these issues have thus been proposed in different fields (biology, image processing, etc.) and although there is an urgent need for changing practices in computer science, how we should proceed is not yet clear. We have recently written an article [3] that contains a completely replicable analysis of the large set of experiments. We took care of doing these experiments in a clean, coherent, wellstructured and reproducible way that allows anyone to inspect how they were performed.

In this presentation, we propose to explain and demonstrate the basics of our workflow and the technology we used.

Context

The reproducible research [3], of which we will illustrate the internals, was done in the context of High-Performance Computing (HPC). Modern HPC infrastructures are typically made of hybrid machines with several CPUs and GPUs. To efficiently exploit all these heterogeneous resources, programmers have to explicitly allocate computation-intensive kernels and manage data transfers between the different processing units. It is common to rely on task-based runtimes that provide the right abstraction to efficiently perform such tedious management and optimizations. Designing and configuring such runtime systems is itself a challenging problem for which we believe the systematic use of modeling and simulation can provide an answer. To support our claim, we have ported StarPU [1], a dynamic task-based scheduling runtime for hybrid architectures, on top of the Sim-Grid simulator [2] and we have shown that this combination allowed to provide accurate performance predictions for dense linear algebra kernels. To validate our approach and our models, we had to perform measurements on several, sometimes not dedicated, machines, and work on complex beta code that often needed to be modified to fit our needs. In this context, even the smallest misunderstanding or inaccuracy about a parameter at small scale often results in a completely different behavior at the macroscopic level. It was thus crucial to carefully collect all the useful metadata, to use a well-planed experiment design and a coherent analysis, all in order to be able to easily reproduce the experiment results.

From Experiments to Article: Reproducibility vs. Replicability

We can distinguish between two main trends in reproducible research. The first one aims at completely automatizing the whole process (experiments, analysis and the final article), allowing others to replicate it. This approach is particularly popular in computational sciences. However, maintaining such process functional in time and portable across platforms is quite burdensome. Additionally, it imposes restrictions on what can be done, it is difficult to extend and hence not always usable on a daily basis. Furthermore, focusing on automatic replication hides information on why and how things were done. The second approach is much lighter and relies on literate programming of code and analysis scripts (i.e. documenting why and how to use the code) and on the systematic use of recipes (i.e., how to build environments). These two approaches apply to both experiments and analysis part of the research.

Reproducible Experimentation

When measuring execution time of a complex application on a modern computer, two runs in a row with the same setup rarely provide the exact same timing in nanoseconds. Such experiments are de facto not replicable, but only reproducible provided that enough details about the configuration are available. We have thus taken a particular care to ensure that all our experiment results are reproducible and will describe in our presentation the most important ones.

Before each measurement, we automatically log information about the machine, such as the current CPU frequency and governor, the memory hierarchy, the versions of Linux and gcc, etc. We also systematically recompile software before using it to conduct experiments, and keep all the configuration and compilation outputs. Additionally, we enforce that all changes to the source code are committed in the revision control system and we keep track of the hash of the Git/SVN version of the source code. This enables provenance tracking, i.e, to know which results were obtained with which code, so they can be later easily investigated, compared or reproduced. Finally, we save the measurement results (makespan, GFlop/s rate) together with execution traces. During the development period, workflow and data format went through several changes and adjustments and we used a laboratory notebook to keep track about all modifications.

Replicable Analysis and Article

Unlike experimentation, analysis is only concerned with the final output (figures, tables, numerical summaries) and not with the time it takes to compute it. Therefore, for a given experimental data set, results of the analysis should always remain the same, which makes them by essence replicable. Reproducible research tools and technologies are still in their infancy and there are many different alternatives (Rstudio/knitr, ActivePapers, Elsevier executable papers, etc.). We tried several of them and we have finally chosen to use org-mode, more precisely orgbabel, that allows to combine multiple processing and analysis codes as well as their output results inside a plain text file. Compared to other alternatives, we feel that org-mode is mature, flexible and lightweight enough for both a daily usage and to allow to produce high quality reproducible articles. The article we have written combines within a single plain text file the body of the article along with all the analysis, all of which are later exported into standard pdf document. It has a hierarchical structure, with different types of code, including Shell (to manipulate data files), R (for plotting figures) and L A T E X (to finely control formatting details). We made all raw data and traces publicly available on figshare [4]. To replicate the article, all these data should be downloaded and unpacked. Then all useful information is extracted into csv files by parsing and filtering raw data files that contain many additional metadata. Finally, R is used to load, process and plot data, ensuring that all figures are consistent throughout the whole article. This way, our experimental and analysis results can be inspected, referenced, or even reused in any other research project.

Conclusion

The approach we propose to present was very effective for doing reproducible research in our context. We strongly believe that it can be easily applied to many other situations and that presenting such methodology can be beneficial to the audience.

Introduction

Les -Définition utilisant les distributions :

X ≤ st Y ⇔ ∀i, 1 ≤ i ≤ n, n k=i d2[k] ≤ n k=i d1[k].
On utilisera de façon équivalente X ≤ st Y et d2 ≤ st d1.

L'ordre ≤ st nous intéresse particulièrement, car il permet de borner les récompenses croissantes, cette propriété est souvent vraie pour les récompenses utilisées en évaluation de performances : moyenne, moment d'ordre quelconque, etc.

Méthode de réduction de la taille d'une distribution en la bornant [1]

Nous

Modèle d'une file FIFO

Soit A(k) une variable aléatoire représentant le trafic entrant dans la file pendant le k-ième intervalle temporel (slot). Nous notons respectivement par Q(k) et D(k) la longueur du tampon et le nombre de sorties pendant ce même slot. Soit B la capacité du tampon et S la capacité de service. On se place dans un modèle où les arrivées ont lieu avant les services. Les services se terminent en fin d'intervalle temporel. Un paquet reste donc au minimum 1 slot dans la file. L'équation de récurrence sur la longueur du tampon dans ce modèle est bien connue [3] :

Q(k) = min(B, (Q(k -1) + A(k) -S) + ), (1)
où (X) + = max(X, 0). On obtient également la distribution de sortie grâce à cette récurrence, pour tout k :

D(k) = min(S, Q(k -1) + A(k)). (2)
L'équation 1 (resp. l'équation 2) définit une chaîne de Markov en temps discret si les arrivées sont indépendantes et le processus A(k) est stationnaire.

Bornes sur la file liées aux bornes sur les histogrammes de trafic

Pour une file à capacité finie B, nous injectons les trafics bornants (inférieurs et supérieurs) et nous établissons un résultat de comparaison stochastique. Supposons que le trafic d'entrée soit Ã(k), un processus d'arrivée qui borne supérieurement ou inférieurement le trafic réel, on note Q(k) l'état de la chaîne de Markov associée à ce processus d'arrivée. On note D(k) le nombre de sorties sous les mêmes hypothèses. Le théorème 4.3.9 de [4] permet d'établir les propriétés suivantes :

Proposition 1 (Bornes supérieures) Si A(k) ≤ st Ã(k), ∀k ≥ 0 alors Q(k) ≤ st Q(k), ∀k ≥ 0 Proposition 2 (Bornes inférieures) Si Ã(k) ≤ st A(k) pour tout k ≥ 0, alors Q(k) ≤ st Q(k) pour tout k ≥ 0.
D'après l'équation 2, nous avons également des bornes sur le processus de sortie.

Proposition 3 Si A(k) ≤ st Ã(k) pour tout k ≥ 0, alors D(k) ≤ st D(k), ∀k ≥ 0. De même, si Ã(k) ≤ st A(k) ∀k ≥ 0, alors, D(k) ≤ st D(k), ∀k ≥ 0.
Et puisque ces propositions sont vraies pour toutes les dates k, elles sont également vraies pour les versions stationnaires qui existent lorsque les chaînes sont ergodiques.

Conclusion

Notre méthode de bornes stochastiques offre un compromis intéressant entre la précision et la vitesse de calcul. En effet, selon le choix de la réduction apporté à la taille des distributions, nous pouvons déterminer des bornes pertinentes et un encadrement très fin du résultat exact en des temps relativement courts. Nous montrons également que l'élément du réseau est stochastiquement monotone ce qui nous permet d'assurer que les distributions bornantes calculées représentent bien des bornes de la distribution exacte. Ainsi, en appliquant nos algorithmes de bornes sur l'histogramme représentant la trace de trafic en entrée, nous définissons des histogrammes bornants sur les différentes mesures de performance : les probabilités de blocage, l'espérance de la longueur du tampon, etc. Definition 1 [2] The stationary sequence {X n } n≥1 , with M n = max{X 1 , . . . , X n } and F as the marginal distribution function, is said to have the extremal index θ ∈ [0, 1] if for each 0 < τ < ∞ there is a sequence of real numbers (thresholds)

u n = u n (τ) such that lim n→∞ n(1 -F(u n )) = τ and (1) lim n→∞ P{M n ≤ u n } = e -τθ
hold.

Some of the interpretations of θ are [2] :

• Relation to the mean cluster size : A cluster is considered to be formed by the exceedances in a block of size r n (r n = o(n)) in n with cluster size ξ = rn i=1 1(X i > u n ) when there is at least one exceedance within r n . The point process of exceedances which counts the number of exceedances, normalized by n, of {X i } n i=1 over threshold {u n } converge weakly to a compound Poisson process (CP) with rate θτ under condition (1) and a mixing condition, and the points of exceedances in CP correspond to the clusters . Then θ = (Eξ) -1 .

• We also have P{M n ≤ x} = F nθ (x) + o(1), n → ∞. Hence θ allows us to evaluate a limit distribution of the maximum of the node degree and it is also helpful to find quantiles of the maxima. These quantiles show the large degrees in the network which arise with a certain probability. The main contributions in this work are as follows.

We study the extremal and clustering properties due to degree correlations in large graphs. Since the network under consideration is known only through Application Programming Interfaces (API), different graph exploring or sampling algorithms are employed to get the samples. The considered algorithms are Random Walk based that are widely discussed in the literature (see [1] and the references therein). In order to facilitate a painless future study of correlations and clusters of degrees in large networks, we propose to abstract the cluster statistics to a single and handy parameter, θ. We derive analytical expressions of θ for different sampling techniques. Finally different estimators of θ are tried out and several other applications are proposed.

Model and Algorithms

We consider networks represented by an undirected graph G with N vertices and M edges. In accordance with the data from most of the real networks, we assume the network is not known completely (with N and M unknown) and also assume correlation in degrees between neighbor nodes. The dependence structure in the graph is described by the joint degree-degree probability density function f(d 1 , d 2 ) (see e.g., [3]). The probability that a randomly chosen edge has the end vertices with degrees

d 1 ≤ d ≤ d 1 + ∆(d 1 ) and d 2 ≤ d ≤ d 2 + ∆(d 2 ) is (2 - δ d 1 d 2 )f(d 1 , d 2 )∆(d 1 )∆(d 2 ). Here δ d 1 d 2 = 1 if d 1 = d 2 , zero otherwise. The degree distribution f d (d 1 ) can be calculated from the marginal of f(d 1 , d 2 ) as f(d 1 ) = d 2 f(d 1 , d 2 ) = d 1 E[D] f d (d 1 )
,

where E[D] denotes the mean node degree. Then E

[D] = (f(d 1 , d 2 )/d 1 ) d(d 1 )d(d 2 ) -1
. Most of the results in this paper are derived assuming continuous probability distributions for f(d 1 , d 2 ) and f d (d 1 ) due to the ease in calculating θ for continuous case. In particular, for analytical tractability and with the sup-port of empirical evidences, we assume the bivariate Pareto model for the joint degree-degree tail function.

Description of random walks

The different graph exploration algorithms considered in the paper are Random Walk based and transition kernels are defined for degree state space unlike in previous works where they were defined and well studied for vertex set (see [1] and the references therein). We use f X to represent the probability density function under the algorithm X .

Standard random walk (RW)

In a standard random walk, the next node to visit is chosen uniformly among the neighbours of the current node. Using "mean field" arguments, the joint density function of the standard random walk is derived as f RW (d t+1 , d t ) = f(d t+1 , d t ).

PageRank (PR)

PageRank is a modification of the random walk which with a fixed probability 1c samples random node with uniform distribution and with a probability c, it follows the standard Random walk transition. Its evolution can be described as

f PR (d t+1 |d t ) = cf RW (d t+1 |d t ) + (1 -c)f d (d t+1 ).
Unfortunately, according to our knowledge, there is no closed form expression for the stationary distribution of PageRank and it is difficult to come up with an easy to handle expression for joint distribution. Therefore, along with other advantages, we consider another modification of the standard random walk.

Random walk with jumps (RWJ)

This algorithm follows a standard Random walk edge with probability d t /(d t + α) and jumps to an arbitrary node uniformly with probability α/(d t + α), where d t is the degree of current node and α ∈ [0, ∞] is a design parameter ( [1]). This modification makes the underlying Markov Chain time reversible, significantly reduces mixing time, improves estimation error and leads to a closed form expression for stationary distribution. The joint density for the random walk with jumps is derived as

f JP (d t+1 , d t ) = E[D]f(d t+1 , d t ) + αf d (d t+1 )f d (d t ) E[D] + α .

Calculation of θ

The extremal index of the random walk based sampling algorithms can be calculated by means of copula density as

θ = lim x→1 x -C(x, x) 1 -x = C ′ (1, 1) -1,
where C(u, u) is the Copula function ([0, 1] 2 → [0, 1]) and C ′ is its derivative [4]. With the help of Sklar's theorem, C(u, u) = P X (D 1 ≤ F -1 X (u), D 2 ≤ F -1 X (u)) where X can be RW, PR or RWJ with F -1 X () is the inverse of the stationary distribution function of the corresponding random walk.

Results

The extremal index is calculated analytically for the considered algorithms. Closed form expressions for RW and RWJ are obtained and a lower bound of θ for PR is derived. As for numerical experiments, a random graph is generated as follows : f(d 1 , d 2 ) is taken as a bivariate Pareto distribution with suitable parameters and the degree distribution is calculated accordingly from (2). Then an uncorrelated graph is generated using the configuration model with this degree distribution. Finally the Metropolis algorithm is applied to rearrange the edges in order to get the desired joint degree-degree distribution. The mean field argument for the transition kernel of the RW is checked with the generated graph and found to yield a reasonable match. θ is estimated for the different algorithms using the rank estimator of Copula function. Finally, for the application of cluster classification, we derive the number of nodes to be sampled to achieve a particular mean number of clusters.

Introduction

In dynamic networks, the topology of the nodes in a network or links connecting these nodes changes with time. This can be due to mobility, link failure, or nodes failure. We are interested in this paper to study consensus on dynamic networks. Most of the work on consensus in dynamic network settings consider fixed number of nodes that are trying to reach agreement in the presence of either mobility or non-robust links (so only the links are dynamic) [1]. Average consensus on these networks is modeled following a random adjacency matrix A n×n (k) where n is the number of nodes and is fixed while the elements of this matrix are random (being for example i.i.d. at every iteration k). The study of consensus on that model is reduced to studying the convergence of the backward product of random matrices. Some papers give sufficient conditions on the weight matrices at every time iteration that guarantee convergence, others use coefficient of ergodicity as a tool to show the convergence of their system to consensus [2]. However, little study has been made on networks with dynamic number of nodes. In the latter case, the dimensions of the adjacency (and weight) matrices can be unbounded, and thus the traditional tools for studying the consensus are not applicable. We refer in this report to this type of networks : nodes arrive and leave as in a queuing system. We would like to study the effect of averaging in these networks, and to see if the nodes could actually reach consensus.

Model

Networks with dynamic nodes are characterized by nodes' arrival and nodes' departure (see Fig. 1). Each node arrives with a random value X k (where k is the number of arrival and each node is labeled by its arrival number). We suppose that there is initially a node having a label 0 in the system that does not leave (it arrives at time t = 0 and does not depart). Let X 0 , X 1 , X 2 , . . . be pairwise independent random variables following some distribution. Let A k be the arrival time of node k and D k its departure time. At any instance t, nodes within the system are connected to each other by some network topology. Let V(t) be the set of all nodes in the system at time instant t, i.e., V(t) = {k | A k ≤ t < D k }. Let Y k (t) be the value of node k that changes with time depending on its history of connection with other nodes, for example it should be clear that Y k (t) has a constant value before arrival (Y k (t) = X k for t < A k ) and after departure (Y k (t) = Y k (D k ) for t ≥ D k ) but this value changes during the time it spends in the system because of interaction with other nodes. The nodes perform a continuous time averaging :

ẋ(t) = -L t x(t), (1)
where x(t) is the state vector of the nodes present in the system at time t, and x i (t) = Y γ(i) (t) where γ(i) is the arrival number of the i-th oldest node in the system (notice that x 1 (t) = Y 0 (t) is always true for the node labeled 0 because it is the oldest node in system and does not depart). L t is the Laplacian of the graph at time t, and ẋ(t) = ∂x(t) ∂t . We call this model a consensus queue due to its similarity to queuing systems. If the inter-arrival and inter-leave times are exponentially distributed random variables with a FIFO discipline (the nodes first to come are the first to leave except for node 0), then the system is M/M/1 consensus queue with arrival rate λ a and departure rate λ d . The Laplacian of the graph L t is used to give a general model for different graph topologies. However, this work is just preliminary and we will give some simplifications in the next section : 1) we only consider two graph topologies, the complete graph and the tree, 2) the averaging is faster than the dynamics of the queue (i.e., equation (1) converges before the arrival or the departure of a new node). The model described here is interesting because it can be applied to different and diverse applications. Queuing consensus can be for example a model for human interactions and their behavior. Consider a system where people arrive at an open market and products' prices are not fixed. Each customer has an initial estimation of the price of the product that varies depending on the interaction process with other customers in the system. You can also consider this model as a representative of a wireless sensor network where sensors monitor some environmental measurement (as temperature or pressure) where nodes can fail according to a Poisson process and new nodes are added to the network. We are interested then by the average in the system, mainly by :

Z(t) = 1 N(t) N(t) i=1 x i (t), (2)
where N(t) = |V(t)| is the number of nodes that are present in the system at time t and x i (t) are their estimates. We note that Z(t) is not a continuous process in general because with every arrival (or departure) with an estimate X k (or Y k ) different from Z(t), the process jumps to a different value.

Simple Network Topologies

There are two sources of randomness in this model, the first one is the input estimate X that follows some distribution, and the other one is the queuing system with random arrivals and departures. In the following sections, we will characterize the average in the system and the output process by considering several simplifications : 1. Complete Graph : the network is a full mesh network (all nodes in the system are connected to each others) and once a node enters a system, all nodes will have the average of nodes presented (instantaneous averaging), i.e., x i (t) = Z(t) for all i in the system. 2. Directed Tree : nodes arriving can only connect to one node chosen uniformly at random in the system, and their estimate changes only once till they leave the network depending on the chosen node's estimate and it's distance from the root.

Complete Graph

Let Z k be the value of Z(t) just after the k-th arrival and before the k+1-th arrival. Then Z k can be written as a weighted average of the nodes in the network, i.e., Z k = k i=0 w i X i , where w i is the weight given to the value of node i and it is a random variable depending on the stochastic arrivals and departures. It is important to study these weights to see how the system preserves the history of old values. To do this we take two extreme cases. The first case of no departures (if λ d = 0), then we have

w i = 1 k + 1
for i = 0, . . . , k.

The second case of very fast departures (λ d >> λ a ), then we suppose each node that arrives, averages with node 0 and then leaves the system, so

w i = ( 1
2 ) k+1-i for i = 1, . . . , k, (

2 ) k for i = 0.

The results are interesting as they show that if the system's departure rate is fast, then the weights for old values decrease exponentially, but if the departure is slow then the weights for old values decrease linearly in k. For future work, we would like to characterize the decrease in the average weight of the history as function of the performance parameter of the queue (as function of ρ = λ a /λ d ).

Directed Tree

For this topology, we assume that there are no departures from the system. Each node arrives, connects uniformly at random to one node (we call it its parent) in the network. Let L be the distance from a newly connected node to the root (node 0). We suppose that each node i arrives at level L, connects to its parent j, and then averages as follows :

Y i (t) = 1 L X i + L -1 L Y j (t) for t > A i .
Notice that Y i (t) = 1 L s∈P X s is a constant for t > A i where P is the set of nodes on the path from the root to i. Since L converges to log n in probability as n → ∞ [3], we conclude that if X 1 , X 2 , . . . are i.i.d random variables of mean µ, then the value Y i (t) for t > A i converges in probability to µ as i → ∞.

c

  i=1 x i ), . . . , (K, M)).La fonction f est une bijection entre S et Π(D).

FIGURE 1 -

 1 FIGURE 1 -Chemins dans un diagramme complet à 5 files et 4 clients avec C 1 = C 5 = 2, C 2 = C 4 = 1 et C 3 = ∞. Le chemin en pointillés représente l'état x = (0, 1, 2, 0, 1).

FIGURE 2 -

 2 FIGURE 2 -T 2,3 (D)

Algorithm 2 :

 2 Simulation parfaite avec D Data: (

Figure 1 )FIGURE 1 -

 11 FIGURE 1 -Plan du port du Havre

FIGURE 2 -

 2 FIGURE 2 -Démarche ECOGRASIM

FIGURE 1 -

 1 FIGURE 1 -Simulation d'une période de 40 itérations sur 4 coeurs. La durée totale de simulation est 25 (temps de couplage entre trajectoires 2 et 4) au lieu de 39 en séquentiel.

  itératons maximales en fonction du nombre de coeurs nombre de coeurs nombre d'itérations

FIGURE 2 -

 2 FIGURE 2 -Modèle de file M/M/1/50 charge ρ = 80%, période de simulation T = 100000. 300 répétitions d'expériences par point. Intervalles de confiance inférieurs à 1%

FIGURE 2 -

 2 FIGURE 2 -Analyse comparative du benchmark NAS LU.A sur les sites de Grenoble et de Lille. Trois compromis distincts sont représentés. Les couleurs représentent la proportion des différents états MPI pour chaque agrégat ((MPI_Recv en bleu, MPI_Wait en rouge et MPI_Send en vert).

References 1 .FIGURE 2 -

 12 FIGURE 2 -Throughput distribution after an intervention on the retransmission score

Figure 1 :

 1 Figure 1: Hit ratio of the original request sequence versus the result of semi-experiment (3).

  ), the curves do not differ considerably. This means that catalog arrivals and individual document request sequences can be modeled as homogeneous Poisson processes. For conciseness, we show in Figure1only the results of semi-experiment (3). Taking into account the previous insights, we propose a two layered model. In the first layer, we consider a Poisson process Γ of rate γ, hereafter called the catalog arrival process. This process models the publications of documents in the catalog. The second layer consists in the document request processes. Specifically, when a document d arrives to the catalog at time a d , we model its request process as an homogeneous Poisson process of rate λ d on the inter-

Figure 2 :

 2 Figure 2: Top: The boxes represent the lifespan and popularity. Bottom: Document request processes. Their superposition is the total request process.

Figure 3 :Proposition 2

 32 Figure 3: Fitting for the Che estimation

Proposition 3

 3 For each setting k ∈ {F, P, N} in the model with N relays, R (k) = C 1 + NC 2 , where C 1 = C d and

FIGURE 1 -

 1 FIGURE 1 -The network model.

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014 Programme Mercredi 11 juin 2014

	12h00-12h30 : Arrivée des participants à l'Inria, espace muséal du bâtiment Euler
	12h30-14h00 : Repas buffet froid, hall du bâtiment Euler
	14h00-15h00 : Exposé de synthèse 1
	Pierre L'Ecuyer
	Challenges in the Stochastic Modeling of Service Systems : Illustrations with Call
	Centers	(p. 5)
	15h00-16h00 : Session 1 : Échantillonnage Parfait/Perfect Sampling
	Stéphane Durand
	A perfect sampling algorithm of random walks with forbidden arcs	(p. 7)
	Anne Bouillard, Ana Bušíc, Christelle Rovetta
	Simulation parfaite dans un réseau fermé de files d'attente	(p. 9)
	16h00-16h30 : Pause café
	de Petri	(p. 11)
	Abderaouf Benghalia, Jaouad Boukachour, Dalila Boudebous
	Contribution à la modélisation et à l'évaluation de la performance du transfert de con-
	teneurs maritimes	(p. 13)
	Marion Dalle, Jean-Marc Vincent, Florence Perronnin
	Catch me if you can !	(p. 15)
	Jeudi 12 juin 2014
	09h00-10h00 :

16h30-18h00 : Session 2 : Simulation Maïder Estécahandy, L. Bordes, S. Collas, C. Paroissin Méthodes d'accélération de Monte-Carlo adaptées à des modèles simples en Réseaux

Session 3 : Analyse de Données/Data Analysis

	Damien Dosimont, Guillaume Huard, Jean-Marc Vincent
	Agrégation temporelle pour l'analyse de traces volumineuses	(p. 17)
	Hadrien Hours, Ernst Biersack, Patrick Loiseau
	A causal study of an emulated network	(p. 19)
	10h00-10h30 : Pause café/viennoiseries
	10h30-11h30 : Exposé de synthèse 2
	Thomas Bonald
	Réseaux de Whittle et applications	(p. 6)
	11h30-12h30 : Session 4 : Systèmes de Caches/Cache Systems
	Felipe Olmos
	The performance of a LRU cache under dynamic catalog traffic	(p. 21)

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014

	Nicaise Choungmo Fofack
	Approximate models for cache analysis with correlated requests	(p. 23)
	12h30-13h00 : Trajet en Pullman jusqu'au restaurant
	13h00-15h00 : Déjeuner à Bijou Plage (Juan-Les-Pins)
	15h00-18h00 : Promenade/découverte de Juan-Les-Pins/tour du Cap d'Antibes
	18h00-: Programme libre
	Vendredi 13 juin 2014
	09h00-10h00 : Session 5 : Contrôle des Réseaux/Network Control
	Tatiana Seregina
	Reward-based Incentive Mechanisms for Delay Tolerant Networks	(p. 25)
	Maialen Larrañaga, U. Ayesta, I.M. Verloop
	Whittle's index in a multi-class queue with abandonments	(p. 27)
	10h00-10h30 : Pause café/viennoiseries
	10h30-11h30 : Exposé de synthèse 3
	Anne Bouillard
	Algorithms for and from network calculus	(p. 6)
	11h30-12h00 : Session 6 : Méthodologie/Methodology
	Luka Stanisic, Arnaud Legrand
	Good Practices for Reproducible Research	(p. 29)
	12h00-13h30 : Repas plateau à la cantine de l'Inria
	13h30-15h00 :

Session 7 : Modèles Stochastiques/Stochastic Models

	Farah Ait Salaht, Jean-Michel Fourneau, Hind Castel, Nihal Pekergin
	Approche par bornes stochastiques et histogrammes pour l'analyse de performance des
	réseaux	(p. 31)
	Konstantin Avrachenkov, Natalia M. Markovich, Jithin Sreedharan
	On Distribution and dependence of extremes in PageRank-type processes	(p. 33)
	Mahmoud El Chamie, Giovanni Neglia, Konstantin Avrachenkov
	Averaging on Dynamic Networks	(p. 35)
	15h00-15h30 :

Table ronde : Logiciels dédiés à la modélisation avec les chaînes de Markov

	15h30-16h00 : Conclusions et fin de l'atelier

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014 Exposés de synthèse

	Pierre L'Ecuyer
	Université de Montréal, chaire internationale Inria
	Challenges

in the Stochastic Modeling of Service Systems :

  Let us consider a random walk in a grid with no forbidden arcs. Let C d be its coalescence time. (i) The number of simulated step before obtaining the coupled point is at most four time the coalescence time. (ii) The expected coalescence time satisfies E[C d ] = O N 2 d log d .

	Lemma 4. If forbidden arcs are chosen randomly, uni-
	formly among all arcs in the grid and if k is the expected
	number of forbidden arcs, then the maximal size |E| of the
	set E, is bounded in expectation : E[ |E| ] ≤ kN π 2 + O(kd +

kN d ).

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014 3. Extension de MCT (EMCT)

	Comme indiqué dans l'introduction, nous proposons
	une extension de MCT qui s'adapte à des cycles de vie
	répétés où l'événement d'intérêt est soit direct, soit en
	conflit avec p autres événements.
	Principe 3.1 Le principe d'EMCT, pour estimer PFD(t),
	est de conditionner chaque événement en conflit de durée
	ξ (1) , . . . , ξ (p) de sorte que la durée de l'événement la plus
	courte soit inférieure à un délai positif ∆ notons ξ(k) , pour k = 1, . . . , p, les durée conditionnées t . Nous
	correspondantes. Leur distribution est donnée, pour tout
	k ∈ {1, . . . , p}, par :
	, X 2 (t), . . . , X m (t)) ∈ {0, 1} m
	est l'état des m composants du système à l'instant
	t ;
	-Φ est une fonction de {0; 1} m dans {0; 1}, type
	fonction de structure, où 0 représente l'ensemble
	des états indésirés du système et 1 les états com-
	plémentaires.
	Nous avons formalisé la PFD pour trois différentes
	structures de composants.
	1. un composant qui tombe en panne selon une
	loi de probabilié quelconque et est réparé du-
	rant des campagnes d'inspections ;
	2. un composant qui a p pannes différentes en
	conflit. Ces pannes apparaissent et sont ré-
	parées selon des lois quelconques. La panne
	d'intérêt est la première panne ;
	3. un composant qui possède une panne de type
	1 en conflit avec une panne de type 2. La panne d'intérêt est la 1 ère panne.

TABLE 1 -

 1 Analyse comparative de la structure 2

	Nb de comp.	1 PFD(8759)=2.9 × 10 -5	2 PFD(8759)=8.2 × 10 -10
	en //	PFD(8759)	Durée	PFD(8759)	Durée
		σ(8759)	(sec)	σ(8759)	(sec)
	MC	2.2 × 10 -5 4.7 × 10 -6	6294	0 0	8489
	EMCT	2.5 × 10 -5 3.7 × 10 -6	6397	0 0	22996
	MC-DISS			8.1 × 10 -10 2.2 × 10 -10	15472
	EMCT-DISS			7.5 × 10 -10 1.6 × 10 -10	25097
	Bibliographie
	1. IEC 62551 ed1.0, "Analysis techniques for depend-
		ability -Petri net techniques", International Elec-
		trotechnical Commission, Geneva, 2012.
	2. B. Tuffin, "Introduction to rare event simulation",
		Inria Rennes (France), AEP9, 2008.
	3. R. Garnier, "Une méthode efficace d'accélération de
		la simulation des réseaux de Petri stochastiques ",
		Thesis, Bordeaux 1 university (France), 1998.
	4. E.E. Lewis, F. Böhm, "Monte Carlo simulation of
		Markov unreliability models", Nuclear Engineering
		and Design, vol. 77, issue 1, pp. 49-62, 1984.

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014 Contribution à la modélisation et à l'évaluation de la performance du transfert de conteneurs maritimes

	Abderaouf Benghalia, Jaouad Boukachour, Dalila
	Boudebous
	Université du Havre
	IUT du Havre, 5 rue Boris Vian, BP 4006,
	76610 Le Havre, France
	abderaoufb@yahoo.fr,
	jaouad.boukachour@univ-lehavre.fr,
	dalila.boudebous@univ-lehavre.fr
	Face au flux croissant de conteneurs et aux con-
	traintes de compétitivité de plus en plus rigoureuses,
	toute plate-forme portuaire doit assurer sa croissance
	et sa rentabilité, tout en sachant maîtriser ses impacts
	environnementaux [1]. A cet effet, le Grand Port Mar-
	itime du Havre (GPMH) a entrepris la construction
	d'un terminal multimodal (cf.

  en raison des coûts associés et de l'impact sur les capacités de manutention de conteneurs[1]. Notre travail concerne la modélisation et la simulation du transfert des conteneurs par navettes ferroviaires entre le futur terminal multimodal et les terminaux maritimes du port du Havre. Il s'agit d'étudier la performance des transferts des conteneurs par rapport au coût, au délai et à la réduction de l'émission de CO 2 . Notre contribution fait état, d'une part, de l'application de la méthode ECOGRAI[2] pour l'identification des indicateurs de performance, et d'autre part, de la simulation à événements discrets pour calculer et mesurer ces indicateurs. A cet effet, nous proposons une démarche dénommée « ECOGRAISIM » [3] qui est basée sur la méthode ECOGRAI et la simulation. La méthode ECOGRAI [2] permet de concevoir et de développer les Systèmes d'Indicateurs de Perfor-

mance (SIP) pour les entreprises industrielles ou de services. Elle permet de guider la conception et l'implantation d'un SIP. Notre démarche ECOGRAISIM consiste à appliquer les quatre premières étapes de la méthode ECOGRAI (cf. Figure

2

). La première étape permet d'établir la grille GRAI (Graphe à Résultats et Activités Inter-reliés)

[3]

pour la modélisation de la structure de pilotage du système et l'identification des centres de décision. GRAI, méthode de modélisation d'entreprises a été crée dans les années 1980

[2]

. Cette méthode, peut être appliquée dans divers domaines. Elle permet

[4]

de détecter les points à améliorer et les points forts du système étudié. GRAI est basée sur une grille caractérisée par des lignes décrivant les niveaux décisionnels (stratégique, tactique et opérationnel), des colonnes décrivant les différentes fonctions du système de pilotage et des flèches représentant le type du flux. L'intersection entre les fonctions et les niveaux de décision nous a permis de déterminer les centres de décision. Dans notre cas d'étude, nous nous sommes focalisés sur un centre de décision appelé « optimiser les processus de manutention et de transfert ». Ce centre est en relation avec le centre de décision appelé « dimensionner et ordonnancer les missions en prenant en compte l'émergence des conteneurs à venir » car un dimensionnement et un ordonnancement optimal influencent directement l'optimisation globale. Notre centre de décision « optimiser les processus de manutention et de transfert » est en interaction avec les informations externes afin de prendre en compte les attentes des clients. Dans la deuxième étape, nous définissons les objectifs concernant les centres de décision et nous identifions les variables de décision dans la troisième étape. En effet, L'objectif de notre centre de décision est la minimisation des différentes pénalités concernant les retards, les coûts et l'émission de CO 2 . Il s'agit également de la minimisation de l'écart entre la date de livraison de conteneurs prévue et la date de livraison de conteneurs au plus tard. Enfin dans la quatrième étape, on obtient les indicateurs de performance. L'exigence d'établir un lien entre les éléments des triplets

[2]

« objectif -variable de décision -indicateur de performance » permet de ne retenir que les indicateurs pertinents dans le cadre de la stratégie globale du terminal multimodal

[2]

. Ils doivent être de type SMART

[4]

(S : Spécifique ; bien décrit et compréhensible, M : Mesurable, A : Atteignable, R : Raisonnable, T : Temporel ; lié avec une échelle de temps.) ou encore appelés indicateurs de performance intelligents. Une fois la liste des indicateurs est établie, la prochaine étape consiste à modéliser le système afin de simuler son comportement.

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014

	Parameter	Definition	Min	Max	Avg	CoV
	bw	minimum bandwidth (Mbps)	1	25	7.1	0.91
	delay	propagation delay (ms)	30	180	86	0.48
	queue1	size of R1 buffer (pkts)	10	400	98	1.1
	queue2	size of R2 buffer (pkts)	10	400	100	0.99
	nlac	Narrow Link Available Capacity (Kbps)	12	3.07e3	630	5
	rwin	C1 advertised receiver window (KB)	74	2.23	290	0.65
	bufferingdelay part of the RTT due to queuing delay (ms)	1	6.76e3	120	2.4
	rtt	Round Trip Time (ms)	84	6.91e3	3.1e2	0.99
	timeouts	number of timeouts (units)	0	682	79	1.5
	retrscore	fraction of retransmitted packets (no unit)	0	0.61	3.7e-3	5.1
	p	fraction of loss events (no unit)	0	0.64	3.8e-3	8.4
	nbbytes	number of bytes sent by the server (MB)	6	150	110	0.21
	tput	throughput (Kbps)	6	1.10e3	280	0.81

TABLE 1 -

 1 Summary of Mininet network emulation experiments dataset

  Il est important de noter que les deux supports des distributions bornantes sont issus de H mais ne sont pas en général égaux. C'est l'algorithme qui fixe le support et la distribution. Les distributions d1 et d2 représentent les bornes optimales sur n états pour la récompense positive croissante r.

	3.	cherchons à réduire la taille des distributions que nous allons employer pour modéliser le trafic et la capacité (débit) de service dans les éléments de réseau. Nous fournissons un encadrement au sens stochastique des résultats. Ces encadrements nous permettent de vérifier que les contraintes de Qualité de Service sont satisfaites. Soit une distribution d, définie par un histogramme à N classes (états), on construit deux distributions bornantes d1 et d2 qui sont définies par des his-togrammes à n < N classes. De plus, d1 et d2 sont les meilleures approximations de d au sens d'une fonc-tion de récompense r fournie par le modélisateur. Nous avons présenté dans [1] un algorithme reposant sur la programmation dynamique pour calculer de telles distributions bornantes optimales. Pour être plus formel, pour une distribution d définie sur H (|H| = N), nous calculons deux distributions d1 et d2 définies respectivement sur H i∈H l r[i]d2[i] est minimale dans l'ensemble des distributions sur n états stochastiquement inférieure à d, i∈H u r[i]d1[i] -i∈H r[i]d[i] est minimale dans l'ensemble des distributions sur n états stochastiquement supérieure à d,

u , H l (|H u | = n, |H l | = n) telles que : 1. d2 ≤ st d ≤ st d1,

2. i∈H r[i]d[i] -

10ème Atelier en Evaluation de Performances Sophia Antipolis, du 11 au 13 juin 2014

Propp, J.G. and Wilson, D.B. -Exact sampling with coupled Markov chains & applications to statistical mechanics. -Random Structures & Algorithms 9 (1996), no. 1-2, 223-252.

Nous avons pris des lois exponentielles uniquement afin de faire des comparaisons avec des valeurs analytiques. Exp(λ) signifie loi exponentielle de taux de défaillance λ.

Lauras, M. -Méthodes de diagnostic et d'évaluation de performance pour la gestion de chaines logistiques, Thèse de doctorat, L'institut National Polytechnique De Toulouse, 2004

. 3. Benghalia, A., Boukachour, J., Boudebous, D. -Évaluation de la performance du trafic des conteneurs maritimes, 9th International Conference on Integrated Design and Production, Tlemcen, Algeria

, 2013. 4. Ducq, Y., Doumeingts, G. -La méthodologie grai et les techniques internet : le projet gvc (grai virtual consulting), 3e Conférence Francophone de MOdélisation et SIMulation, Conception, Analyse et Gestion des Systèmes Industriels, MOSIM'01,

-Troyes, 2001. 5. Beaudry, M. -Indicateurs de performance et tableau de bord. Courses, http ://video.coursgratuits.net/5/p-strategieindicateurs-de-performance.php

A. Müller and D. Stoyan. Comparison Methods for Stochastic Models and Risks. Wiley, New York, NY, 2002.

Sophia -du GDR ASR -du GDR RO -du GDR IM -de l'ANR MARMOTE (projet 12-MONU-0019