L’objectif 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.

Simulation de la file LIFO

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)
         )
}

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
    
    #### A few useful local functions 
    run_task = function() { # runs the last task of the waiting list
      if(length(waiting)>0) {
        running <<- waiting[length(waiting)]
        remaining[running] <<- switch(policy,
                                      npmtn = S[running],
                                      pmtn = min(S[running],remaining[running],na.rm=T),
                                      pmtn_restart = S[running],
                                      pmtn_reset = Service(1,typeservice,x,y)
                                      )
        waiting <<- waiting[-length(waiting)]
      }
    }

    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)}
        running <<- NA
      }
      waiting <<- c(waiting,next_arrival)
      next_arrival <<- next_arrival+1 
      if(is.na(running)) { run_task() }
    }

    #### Main simulation loop
    while(TRUE) { 
      # Look for next event
      dt = NA
      if(next_arrival <=n) { dt = min(dt,(t1[next_arrival]-t), na.rm=T) }
      if(!is.na(running))  { dt = min(dt,remaining[running], na.rm=T)   }
      if(is.na(dt)) { break }
      
      # Update state
      t=t+dt
      if(!is.na(running)) {
        remaining[running] = remaining[running] - dt
        if(remaining[running]<=0) {
          t2[running] = t
          running = NA
          run_task()
        }
      }
      if((next_arrival<=n) & (t==t1[next_arrival])) {
        push_task()
      }
    }
    
    t2-t1
}

Bon, vérifions que ça marche.

FileLIFO(n=40,lambda=.3,"exp",1,0,"npmtn")

##  [1] 0.05148 3.02959 0.61412 1.18420 0.55615 0.15448 0.24710 2.01548
##  [9] 3.72483 2.87035 0.85260 1.71917 1.74751 1.42148 1.21966 0.66956
## [17] 1.83934 3.02835 0.33077 0.38702 1.23669 1.93435 0.05708 0.40590
## [25] 0.17965 0.22894 0.40124 2.30411 1.88015 2.01172 3.60713 0.61554
## [33] 3.27207 1.58602 1.58223 1.51957 0.07370 0.29674 1.50641 2.12806

J’avais les idées claires sur ce que je voulais faire et ce qui m’a pris le plus de temps, c’est de réaliser que j’avais besoin de faire des fonctions à effets de bords et comment faire, ainsi que réaliser le fait qu’avoir une variable qui s’appelle T est une très mauvaise idée… J’ai lancé quelques fois à la main sur des petites instances avec des prints pour vérifier que c’était cohérent, mais je dois dire qu’une fois la npmtn bien codée, j’ai écrit les trois autres en moins de 10 minutes et j’avais assez confiance (je n’ai absolument pas eu à revenir dessus).

Étude des temps de service

stypes = rbind(
  data.frame(typeservice="det",x=0,y=0),
  data.frame(typeservice="uni",x=0,y=2),
  data.frame(typeservice="uni",x=.5,y=1.5),
  data.frame(typeservice="exp",x=1.0,y=0),
  data.frame(typeservice="gamma",x=4,y=1/4),
  data.frame(typeservice="gamma",x=.2,y=1/.2))
stypes$label = as.factor(paste(stypes$typeservice,"(",stypes$x,",",stypes$y,")",sep=""))
df = data.frame();
for (stype in stypes$label) {
  d = stypes[stypes$label==stype,]
  df=rbind(df,data.frame(law=stype, 
                         service=Service(40000,typeservice = as.character(d$typeservice), x=d$x, y=d$y)))
}

J’ai choisi 40000 échantillons. C’est un peu beaucoup, mais ça va vite et ça permet d’avoir des histogrammes bien réguliers et des estimations plutôt propres.

dfsummary = ddply(df,c("law"),summarize, mean=mean(service),sd=sd(service),
                  ci = 2*sd(service)/sqrt(length(service)))
arrange(dfsummary, sd)

##             law   mean     sd       ci
## 1      det(0,0) 1.0000 0.0000 0.000000
## 2  uni(0.5,1.5) 0.9981 0.2887 0.002887
## 3 gamma(4,0.25) 0.9962 0.5002 0.005002
## 4      uni(0,2) 0.9972 0.5778 0.005778
## 5      exp(1,0) 1.0064 0.9992 0.009992
## 6  gamma(0.2,5) 0.9863 2.2033 0.022033

Tout va bien, la moyenne empirique est bien de 1. La variance pour la gamma(0.2,5) est vraiment très grande par rapport aux autres. La gamma(4,0.25) a une variance plus faible que l’exponentielle. Elle fait donc un bon “intermédiaire” entre les uniformes et l’exponentielle.

Regardons les les unes à coté des autres.

ggplot(data=df) + geom_histogram(aes(x=service)) + facet_wrap(~law,) + geom_vline(data=dfsummary,aes(xintercept=mean),color="red")

## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

plot of chunk unnamed-chunk-5

Mmmh, on a des échelles très différentes. Entre la déterministe qui vaut toujours 1 et la gamma(0,2.5) qui a une variance énorme… Regardons leur forme en faisant bien attention au fait que l’échelle est propre à chaque histogramme cette fois-ci

ggplot(data=df) + geom_histogram(aes(x=service)) + facet_wrap(~law,scales="free") + geom_vline(data=dfsummary,aes(xintercept=mean),color="red")

## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

plot of chunk unnamed-chunk-6

Pour gamma(0,2.5), la très grande majorité des valeurs sont en dessous de la moyenne mais c’est compensé par des valeurs “extrèmes” relativement rares (l’échelle en x va au delà de 40…). On a donc là, des lois avec des variances assez différentes les unes des autres, ce qui devrait permettre de voir comment la variance impacte les politiques des files.

Étude détaillée de la fille M/M/1

create <- function(n=10000,lambdas = c(0.2,0.4,0.6,0.8),typeservice="exp",x=1,y=NA,
                   policies = c("npmtn","pmtn","pmtn_restart","pmtn_reset")) {
  d <- data.frame()
  for (policy in policies) {
    for (lambda in lambdas) {
#      print(paste(policy,lambda,sep=" ")) # ça, c'est juste pour suivre l'avancement
      r = FileLIFO(n=n, lambda=lambda, typeservice=typeservice, x=x, y=y, policy)
      d = rbind(d,data.frame(lambda=lambda, policy=policy, resp_m = mean(r), 
                             resp_ci = 2*sd(r)/sqrt(length(r))))
    }
  }
  
  d$input_type = typeservice
  d$input_p1 = x
  d$input_p2 = y
  d$label = as.factor(paste(typeservice,"(",x,",",y,")",sep=""))
  d
}

Essayons avec un échantillonnage relativement grossier et regardons si le temps que ça prend est raisonnable.

system.time(df <- create())

##    user  system elapsed 
##   5.352   0.000   5.358

Moins de 10 secondes sur ma machine. Ça me va. Alors, comment se comparent ces différentes politiques de service ?

p = ggplot(df,aes(x=lambda, y=resp_m, color=policy, shape=policy)) + 
  geom_line() + 
  geom_point() + 
  geom_errorbar(aes(x=lambda, ymin=resp_m-resp_ci, ymax=resp_m+resp_ci),
                width = .02) +
  scale_color_brewer(palette="Dark2")+
  labs(title= "Temps de réponse moyen en fonction de la charge en entrée") + 
  xlab("Débit normalisé d'entrée (clients/unité de temps)\n(Temps de service = 1 unité de temps)") +
  ylab("Unités de temps") +  xlim(0,1) +
  geom_vline(xintercept = 1) + geom_hline(yintercept = 1) + theme_bw()
p

plot of chunk unnamed-chunk-9

Wow, visiblement, pmtn_restart explose complètement. Limitons un peu l’échelle des ordonnées pour y voir plus clair:

p + ylim(0, 6)

## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).

plot of chunk unnamed-chunk-10

Tiens, on ne voit aucune différence notable entre npmtn, pmtn, et pmtn_reset… Les intervalles de confiance sont un peu large ceci dit et peut-être qu’on gagnerait à faire des trajectoires un peu plus longues. Quoi qu’il en soit, si on y réfléchit un peu, il me semble que tout ça s’explique. Avec la politique pmtn_restart, quand les tâches se font interrompre, on recommence au début et on a gaché du temps de calcul. Donc, quand on augmente le débit d’arrivée, il y a beaucoup plus de préemptions et donc plus de travail à faire par tâche. Le système sature donc bien plus tôt. Après, pour arriver à déterminer quand exactement… ça parait super dur car le système étant très instable, il est d’autant plus difficile de déterminer son comportement.

Maintenant, en ce qui concerne pmtn et pmtn_reset, on ne voit probablement pas de différence à cause de la loi du temps de service exponentielle. Comme cette loi est sans mémoire, qu’on reprenne où on s’est arrêté ou qu’on retire un temps de service, c’est pareil… Reste la question de pourquoi on ne voit pas de différence entre npmtn et pmtn. Vu que ces deux politiques sont équivalentes du point de vue du nombre de tâches dans le système, elles ont le même temps de service moyen (Little). Du coup, je ne pense pas qu’il y ait non plus de différence avec la FIFO pour ces politiques là.

Étude de la file M/GI/1

Simulons pour différents types de temps de service

system.time(df <- rbind(df, create(typeservice="det"),
                  create(typeservice="uni",x=0,y=2),
                  create(typeservice="uni",x=.5,y=1.5),
                  create(typeservice="gamma",x=4,y=1/4),
                  create(typeservice="gamma",x=.2,y=1/.2)))

##    user  system elapsed 
##  30.364   0.016  30.409

Bon, c’est plus long, mais ça reste tout à fait raisonnable, hein ? Regardons ce que ça donne…

p = ggplot(df,aes(x=lambda, y=resp_m, color=policy, shape=policy)) + 
  geom_line() + 
  geom_point() + 
  geom_errorbar(aes(x=lambda, ymin=resp_m-resp_ci, ymax=resp_m+resp_ci),
                width = .02) +
  scale_color_brewer(palette="Dark2")+
  labs(title= "Temps de réponse moyen en fonction de la charge en entrée") + 
  xlab("Débit normalisé d'entrée (clients/unité de temps)\n(Temps de service = 1 unité de temps)") +
  ylab("Unités de temps") +  xlim(0,1) + ylim(0, 6) +
  geom_vline(xintercept = 1) + geom_hline(yintercept = 1) + theme_bw() +
  facet_wrap(~label)
p

## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 5 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).

plot of chunk unnamed-chunk-12

Alors:

Exp, c’est le même qu’avant, tout va bien.
Det: pmtn_restart et pmtn_reset ont exactement le même comportement dans le cas déterministe (logique…) et saturent bien avant 1 (, ce qui est également logique puisqu’elles font strictement plus de travail que les autres). Par contre, je suis surpris de voir une différence entre pmtn et npmtn (du coup, je ne suis pas sûr que mon argument précédent tienne la route). Oh, mais oui, finalement, ça s’explique très bien. Si je prend juste deux tâches de durée 1 et que leur exécution est suceptible de se chevaucher, le temps de réponse total est bien meilleur avec npmtn qu’avec pmtn car pmtn fait tourner une tâche avec un remaining bien supérieur à celui de la tâche en cours. Le fait que ces deux politiques soient équivalentes dans le cas exp est donc bien lié au côté sans mémoire mais ça ne me parait pas évident à expliquer.
Uni(0,2), Uni(0.5,1.5) et gamma(4,0.25): si on regarde bien chacune des courbes, elles sont bien entre det et Exp. On a l’impression de passer de l’une à l’autre quand la variance augmente. pmtn_restart sature avant tout le monde et pour les autres, on a toujours npmtn \(\leq\) pmtn \(\leq\) pmtn_reset. Et plus la variance du temps de service est faible, plus on voit de différence entre les trois politiques, pmtn restant bien au milieu et semblant peu affectée.
Gamma(0.2,5): alors ça, c’est bizarre. Bon, pmtn_restart explose avant tout le monde et vu la variance énorme qu’il y a, ça n’est pas surprenant. Il suffit de se taper une fois une tâche avec un temps de traitement un peu long pour qu’elle reste indéfiniment dans le système, se faisant à chaque fois préempter et recommencer depuis le début. Par contre, ce qui est amusant, c’ets le comportement de pmtn_reset qui donne des temps de réponse moyens inférieurs à 1! Intuitivement, dès qu’une tâche a un temps de réponse un peu grand, elle se fait préeempter et recommencera avec un temps de réponse super petit. Impossible de garder un grand temps de service, seuls les petits peuvent passer donc le système sature a priori plus loin. Essayons de voir ce que ça donne quand on augmente la charge…

system.time(df <- rbind(df, create(lambdas = seq(from=1,to=5,by=.2),
                                   typeservice="gamma",x=.2,y=1/.2,policies=c("pmtn_reset"))))

##    user  system elapsed 
##   9.400   0.004   9.410

Et regardons ce que ça donne…

p = ggplot(df,aes(x=lambda, y=resp_m, color=policy, shape=policy)) + 
  geom_line() + 
  geom_point() + 
  geom_errorbar(aes(x=lambda, ymin=resp_m-resp_ci, ymax=resp_m+resp_ci),
                width = .02) +
  scale_color_brewer(palette="Dark2")+
  labs(title= "Temps de réponse moyen en fonction de la charge en entrée") + 
  xlab("Débit normalisé d'entrée (clients/unité de temps)\n(Temps de service = 1 unité de temps)") +
  ylab("Unités de temps") +  xlim(0,max(1,df$lambda)) + ylim(0, 6) +
  geom_vline(xintercept = 1) + geom_hline(yintercept = 1) + theme_bw() +
  facet_wrap(~label)
p

## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 3 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 5 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_path).

plot of chunk unnamed-chunk-14

Bon, c’est un comportement pour le moins suprenant. Ça fini par monter un peu (avec des intervalles de confiance assez élevés d’où l’instabilité de l’estimation) mais il faut vraiment une charge très élevée. C’est une file hyper stable…

Éval de perf, DM: Les derniers seront les premiers (RICM4, 2014)

Arnaud Legrand

26/02/2015

Simulation de la file LIFO

Étude des temps de service

Étude détaillée de la fille M/M/1

Étude de la file M/GI/1