Sources : Pour ce DM nous nous sommes inspirés de cette page http://rpubs.com/alegrand/13532 L’objectif de ce DM est d’analyser l’importance de la distribution du temps de service sur le temps de réponse dans une file d’attente M/GI/1 avec un ordonnancement LIFO. Le processus d’arrivée est un processus de Poisson de taux (débit) Les clients ont un temps de service de moyenne 1 pris comme unité de temps de référence.

Pour que l’on s’accorde sur les mêmes termes : Le temps de réponse : est le temps entre l’arrivée et le départ d’un processus (noté R) Le temps d’attente : est le temps durant lequel le processus est préempté (noté A) Le temps de service : est le temps durant lequel le processus s’exécute (noté S). Donc dans le cas d’une politique restart, ce temps peut être plus grand que le temps généré aléatoirement Ainsi on aurait R = A + S

Simulation de la file LIFO :

options(warn=-1)
set.seed(10)
library(plyr)
library(ggplot2)

Service <- function(n=1,typeservice,x,y) {
# genere un temps de service
  switch(typeservice,
         det = rep(1,n),
         uni = runif(n,x,y),
         gamma = rgamma(n,shape=x,scale=y),
         exp = rexp(n,x)
         )
}

Question 1

Nature des lois de service

Nature des lois de service : Illustrons les différences de natures entre les différentes lois de temps de service :

set.seed(10)

laws        = data.frame()
laws        = rbind(laws, data.frame(name="det(0,0)", fun="det", x=0, y=0))
laws        = rbind(laws, data.frame(name="uni(0,2)", fun="uni", x=0, y=2))
laws        = rbind(laws, data.frame(name="exp(1,0)", fun="exp", x=1, y=0))
laws        = rbind(laws, data.frame(name="gamma(4,0.25)", fun="gamma", x=4, y=0.25))
laws        = rbind(laws, data.frame(name="gamma(0.2,5)", fun="gamma", x=0.2, y=5))

n = 100
df = data.frame()

for (law in 1:nrow(laws)) {
  tmp = data.frame(value=Service(n, typeservice=as.character(laws[law,]$fun), laws[law,]$x, laws[law,]$y))
  tmp$id  = 1:length(tmp$value)  
  tmp$name = laws[law,]$name
  
  df = rbind(df, tmp)
}
ggplot(data=df, aes(x=id, y=value)) + 
  geom_bar(data=df,stat="identity", position="identity") + 
  #geom_hline(yintercept = mean(df$value), color="red") +
  facet_wrap(~name) +
  ylab("Time") + xlab("Client") + 
  ggtitle("Temps de service")

plot of chunk unnamed-chunk-3

Afin de se donner un ordre de comparaison entre les différentes lois, on génère côte à côte les temps de service. Une première remarque est la grandeur d’un certain nombre (limité certes) de temps de service de “gamma(0.2,5)” Intéressons-nous de manière plus précise à chacune des lois :

fun = df[df$name=="det(0,0)",]
ggplot(data=fun, aes(x=id, y=value)) + 
  geom_bar(data=fun,stat="identity", position="identity") + 
  geom_hline(yintercept = mean(fun$value), color="red") +
  ylab("Time") + xlab("Client") + 
  ggtitle(paste("Temps de service", as.character(fun[1,]$name), sep=" : "))

plot of chunk unnamed-chunk-4

fun = df[df$name=="uni(0,2)",]
ggplot(data=fun, aes(x=id, y=value)) + 
  geom_bar(data=fun,stat="identity", position="identity") + 
  geom_hline(yintercept = mean(fun$value), color="red") +
  ylab("Time") + xlab("Client") + 
  ggtitle(paste("Temps de service", as.character(fun[1,]$name), sep=" : "))

plot of chunk unnamed-chunk-5

Les lois déterministes et uniformes sont triviales. Leur moyenne de temps de service est de 1 (~0,90 pour la loi uniforme par manque d’échantillons). On peut ajouter que dans le cas uniforme les temps sont compris entre 0 et 2 secondes par choix.

fun = df[df$name=="exp(1,0)",]
ggplot(data=fun, aes(x=id, y=value)) + 
  geom_bar(data=fun,stat="identity", position="identity") + 
  geom_hline(yintercept = mean(fun$value), color="red") +
  ylab("Time") + xlab("Client") + 
  ggtitle(paste("Temps de service", as.character(fun[1,]$name), sep=" : "))

plot of chunk unnamed-chunk-6

Loi exponentielle :

La plupart des temps de service sont très courts (voire proches de zéro surtout quand x grandit) et d’autres (un nombre limité) assez longs : c’est typique de la distribution exponentielle.

fun = df[df$name=="gamma(4,0.25)",]
ggplot(data=fun, aes(x=id, y=value)) + 
  geom_bar(data=fun,stat="identity", position="identity") + 
  geom_hline(yintercept = mean(fun$value), color="red") +
  ylab("Time") + xlab("Client") + 
  ggtitle(paste("Temps de service", as.character(fun[1,]$name), sep=" : "))

plot of chunk unnamed-chunk-7

fun = df[df$name=="gamma(0.2,5)",]
ggplot(data=fun, aes(x=id, y=value)) + 
  geom_bar(data=fun,stat="identity", position="identity") + 
  geom_hline(yintercept = mean(fun$value), color="red") +
  ylab("Time") + xlab("Client") + 
  ggtitle(paste("Temps de service", as.character(fun[1,]$name), sep=" : "))

plot of chunk unnamed-chunk-8

Loi gamma :

Nous avons choisi deux coefficients x={4, 0.25} afin de nous éloigner de la loi exponentielle (avec x=1).

Dans le cas “gamma(4,0.25)”

var(df[df$name=="gamma(4,0.25)",]$value)
## [1] 0.3686
mean(df[df$name=="gamma(4,0.25)",]$value)
## [1] 1.118

On constate que la moyenne des temps de service reste à 1 et que les valeurs tournent autour de 0,2 et 3 secondes.

Dans le cas “gamma(0.2,5)” :

var(df[df$name=="gamma(0.2,5)",]$value)
## [1] 6.561
mean(df[df$name=="gamma(0.2,5)",]$value)
## [1] 1.092

Des temps de service très variables entre 0 et 13 seconde. En effet malgré la moyenne de 1, la variance reflète un écart important des valeurs.

Conclusion

Dans l’ensemble, la moyenne des temps de service générés est autour de 1. Seule leur variance diffère. Ainsi cela nous permettra de tester le comportement et la robustesse de nos politiques.

Voici un schéma résumant la moyenne des temps de service (en barre) mais aussi leur variabilité.

diff = ddply(df, c("name"), summarise, mean = mean(value), variance = var(value))
ggplot(diff, aes(x=name)) + 
  geom_bar(stat="identity", position="identity", aes(y=mean)) +
  geom_point(aes(y=variance), size=5, color="red") +
  geom_line(aes(y=variance))
## geom_path: Each group consist of only one observation. Do you need to adjust the group aesthetic?

plot of chunk unnamed-chunk-11

Question 2

Etude détaillée de la file M/M/1 - LIFO

FileLIFO <- function(n,lambda,typeservice,x,y,policy) {
    # simulates a M/GI/1 LIFO queue with different preemption policy
    # parameters:
    #    n :  total number of jobs
    #    lambda : arrival rate
    #    typeservice : service law (det uni gamma exp)
    #    x ,y : parameters of the service law
    #    policy: npmtn, pmtn, pmtn_restart, pmtn_reset
    # return value:
    #    vector with response time of each task assuming the queue is initially empty
    
    A <- rexp(n,lambda)         # inter arrival
    t1 <- cumsum(A)             # arrival dates
    t2 <- rep(NA,n)             # completion dates
    S <- Service(n,typeservice,x,y) # initial service times
    
    #### Variables that define the state of the queue
    t = 0               # current time
    remaining = rep(NA,n)  # how much work remains to do for each task
    running = NA        # index of the currently running task
    waiting = c()       # stack with tasks which have arrived and have not been completed yet
    next_arrival = 1    # index of the next task to arrive
    sys_vide = 0        # by curiosity let's count the number of times the system is empty
    sys_interrupted = 0 # by curiosity let's count the number of times the system is interrupted
    used = rep(0, n)    # by curiosity let's count the total of time used by a process (including its restart)
    #interrupted = data.frame()  # by curiosity let's save the preemptions for counting them later
    
    #### A few useful local functions 
    run_task = function() { # runs the last task of the waiting list
      #if there is still a process in the stack to run, we take the last one in
      if(length(waiting)>0) {
        running <<- waiting[length(waiting)] #it unstacks th e last one
        #we update the remaning time : it depends on the policy chosen
        remaining[running] <<- switch(policy,
                                      #the process restart from the beginning
                                      npmtn = S[running],
                                      
                                      #At first loop it takes S[running], but after remaining[running] is returned
                                      pmtn = min(S[running],remaining[running],na.rm=T),
                                      
                                      #the process restart from the beginning
                                      pmtn_restart = S[running],
                                      
                                      #the process restart with a shorter time
                                      pmtn_reset = Service(1,typeservice,x,y)
                                      )
        waiting <<- waiting[-length(waiting)] #renvoie le tableau privé de l'element a la position length(waiting) => le dernier element empilé
      }
    }

    push_task = function() { # insert the next_arrival-th task to the waiting list
                             # and run it if there is preemption
      if(policy != "npmtn") {
        if(!is.na(running)) {
          waiting <<- c(waiting,running) #it stacks the current process, if t==t1[next_arrival]
          #interrupted = rbind(interrupted, data.frame(id=running, when=t, remaining=remaining[running]))
        }
        running <<- NA
      }
      waiting <<- c(waiting,next_arrival)
      next_arrival <<- next_arrival+1 
      if(is.na(running)) { run_task() }
    }

    #### Main simulation loop
    while(TRUE) { 
      #dt is the time we will minus the remaning running process time
      # Look for next event
      dt = NA
      #have we reached the number max of process arrival ?
      if(next_arrival <=n) { dt = min(dt ,(t1[next_arrival]-t), na.rm=T ) } #dt=t1[..]-t it's a bit akward to take the min with dt because NA are deleted. Notice that t1[next_arrival]-t=dt is the next process arrival time - current time. So this is the elapsed time
      
      #is some process running ? if so, we finish executing it if its remainning time is lower than the delta time of the next arrival
      if(!is.na(running))  { dt = min(dt,remaining[running], na.rm=T) }
      
      #we've reached the max of process (no next arrival) and no process is running, we take off
      if(is.na(dt)) { break }
      
      # Update state
      t=t+dt
      if(!is.na(running)) {
        remaining[running] = remaining[running] - dt
        used[running] = used[running] + dt
        if(remaining[running]<=0) {
          t2[running] = t
          running = NA
          
          #it unstacks the last process and it launchs it depending on the policy chosen (resume, restart, restart with lower time), and affect remain[running] to the corresponding time service
          run_task()
        }
      }
      if((next_arrival<=n) & (t==t1[next_arrival])) {
        push_task() #insert the next_arrival-th task to the waiting list and run it if there is preemption
      }
      
    }
    
    #count the number of times the system is empty
    for(i in 2:n) {
      if ( t1[i]>t2[i-1]) {
        sys_vide = sys_vide + 1
      }
    }
    
    list (jobs = data.frame(arrival = t1, completion=t2, service=used, theoricService=S, response=(t2-t1)), 
          sys_vide = sys_vide)
}    

Etudions de manière plus précise le rôle de X dans la loi exponentielle :

set.seed(10)
lambda_min  = 0.1
lambda_max  = 0.95
step        = 0.05
lambdas     = c(.2, .4, .6, .8) #seq(lambda_min, lambda_max, step)
n           = 5000
policies    = c("npmtn","pmtn", "pmtn_restart", "pmtn_reset")

laws        = data.frame()
laws        = rbind(laws, data.frame(name="exp(1,0)", fun="exp", x=1, y=0))
laws        = rbind(laws, data.frame(name="exp(5,0)", fun="exp", x=3, y=0))
laws        = rbind(laws, data.frame(name="exp(0.1 ,0)", fun="exp", x=0.2, y=0))

df = data.frame()
for(law in 1:nrow(laws)) {
  for(policy in policies){ 
    for(lambda in lambdas) {
        tmp=FileLIFO(n, lambda, typeservice=as.character(laws[law,]$fun), laws[law,]$x, laws[law,]$y, policy)$jobs
        tmp$mode = policy
        tmp$lambda = lambda
        tmp$id=1:length(tmp$arrival)
        tmp$law = laws[law,]$name
        df = rbind(df, tmp)
    }
  }
}

library(plyr)
res = ddply(df, c("law", "mode", "lambda"), summarize, n = length(response), responseM=mean(response), sd_response=sd(response))
ggplot(data=res, aes(x = lambda, y = responseM, color = mode, shape = mode)) + 
  geom_line() +
  geom_point() + 
  geom_errorbar(width = 0.02, aes(x = lambda, y = responseM, ymin = responseM - 2 * sd_response/sqrt(n), ymax = responseM + 2 * sd_response/sqrt(n))) + 
  geom_vline(xintercept = 1) +
  geom_hline(yintercept = 1) +
  theme_bw() +
  facet_wrap(~law) +
  xlab("lambda") +
  ylab("Temps de réponse moyen") +
  ggtitle("Influence du taux d'arrivée (lambda) sur le temps de réponse")

plot of chunk unnamed-chunk-14

Nous avons testé différents paramètres pour la loi exponentielle générant le temps de service afin de trouver le plus cohérent. Nous pouvons par exemple constater qu’avec un paramètre x trop faible, les temps de service sont trop élevés et les résultats sont difficilement exploitables, tandis que s’ils sont trop faibles nous ne pouvons pas déterminer les caractéristiques de chaque politique puisque les requêtes se font rarement en concurrence.

Fixons nous une valeur X = 1

set.seed(10)
lambda_min  = 0.1
lambda_max  = 0.95
lambdas     = c(.2, .4, .6, .8) #seq(lambda_min, lambda_max, step)
step        = 0.05
n           = 10000
x           = 1
y           = 1
policies    = c("npmtn","pmtn", "pmtn_restart", "pmtn_reset")

df = data.frame()
for(policy in policies){ 
  for(lambda in lambdas) {
      route=FileLIFO(n, lambda, typeservice="exp", x, y, policy)
      tmp=route$jobs
      tmp$mode = policy
      tmp$lambda = lambda
      tmp$id=1:length(tmp$arrival)
      tmp$n = n
      tmp$sys_vide = route$sys_vide
      tmp$sys_interrupted=(tmp$n-tmp$sys_vide)
      df = rbind(df, tmp)
  }
}


library(plyr)
res = ddply(df, c("mode", "lambda", "n", "sys_vide", "sys_interrupted"), summarize, serviceM= mean(service), sd_service=sd(service), responseM=mean(response), sd_response=sd(response))
ggplot(data=res, aes(x = lambda, y = responseM, color = mode, shape = mode)) + 
  geom_line() +
  geom_point() + 
  geom_errorbar(width = 0.02, aes(x = lambda, y = responseM, ymin = responseM - 2 * sd_response/sqrt(n), ymax = responseM + 2 * sd_response/sqrt(n))) + 
  geom_vline(xintercept = 1) +
  geom_hline(yintercept = 1) +
  theme_bw() +
  xlab("lamba") +
  ylab("Temps de réponse moyen") +
  ggtitle("Influence du taux d'arrivée (lambda) sur le temps de réponse")

plot of chunk unnamed-chunk-16

Nous poussons l’échantillonnage à 10 000 arrivées afin d’obtenir une courbe lissée et de conserver des intervalles de confiance petits avec un degré de confiance alpha=95% soit un phi(95%) = qnorm(0,95 - (1-0,95)/2) ~ 2.

Nous pouvons d’ores et déjà constater que les politiques (npmtn), resume (pmtn), et reset (avec un autre temps de service : pmnt) sont stables contrairement à pmtn_restart dont le temps de réponse explose quand le débit d’inter arrivée tend vers 1.

Les intervalles de confiance relativement petits contiennent les temps estimés. Cela atteste de la confiance que l’on peut accorder à cette estimation des temps de réponses.

Dans cette question, le débit d’inter arrivée est gouverné par une loi exponentielle (M/M/1) dont nous avons pu observer par ailleurs à la question précédente que plus lambda est grand, plus sa courbe de densité est aplatie contre les axes x, y donc les valeurs aléatoires tirées sont petites. Des temps d’inter arrivée plus petits génèrent davantage de chevauchements entre processus. Cela augmente le temps de réponse en moyenne.

Avoir la courbe pmtn_restart qui explose est prévisible, puisque chaque processus préempté recommence son exécution à partir de zéro.

Si l’on observe le nombre de fois qu’un processus est préempté en moyenne quand lambda augmente, il n’est guère étonnant que la politique “restart” soit instable

#mean = ddply(res, c("lambda"), summarize, sys_int= mean(sys_interrupted))
#geom_smooth(data=mean, (aes(y=sys_int)), method="lm")
res = res[res$mode != "npmtn",] #la maniere dont on compte le nombre de fois que le système est vide, font que les valeurs pour ce mode sont non significatives. En effet dans un diagramme d'exectution toute les requêtes se font à la suite.
ggplot(data=res, aes(x=lambda) )+
  geom_point(aes (y=(sys_interrupted), colour=(sys_interrupted), shape=mode)) +
  geom_smooth(data=res, (aes(y=sys_interrupted)), method="lm")+
  xlab("Lambda") +
  ylab("Le nombre de préemptions") +
  ggtitle("Le nombre de préemptions selon le débit d'inter arrivée (lambda) ")

plot of chunk unnamed-chunk-17

Par ce graphe nous montrons la relation entre le débit d’inter arrivée et le nombre de préemptions qui dépend indirectement du temps de service. Nous observons ainsi une régression linéaire qui indique une proportionnalité entre ces deux variables.

Question 3

Étude de la file M/GI/1 − LIFO

Traçons le temps moyen de réponse en fonction de