Dans ce DM, nous nous intéressons à la simulation d’une file d’attente de type M/GI/1 avec différents temps de service.

Le système possède différents paramètres :

Simulateur

On reprend le code sur la page du sujet en modifiant les valeurs de retour pour ajouter différentes informations sur la simulation. (pour avoir plus d’information si nécessaire par la suite).

La valeur de retour sera une data frame, qui pour chaque entrée aura les champs suivants :

Note : Nous nous sommes rendus compte que finalement, nous n’avions pas besoin de tous les champs retournés, nous les avons donc commenté pour gagner en performance. De plus, nous réalisons une modification dans le simulateur pour faire tourner à vide les premiers essais (phase d’initialisation du système) et ne garder que les derniers pour avoir une meilleure variance.

set.seed(42)
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
    #     data frame with values :
    #       - inter : interarrival time
    #       - arr : arrival date
    #       - dep : leave time
    #       - serv : service time
    #       - tps : time spent in the system
    #       - att : tps-serv (ideally 0 for each value)
  
    ####### HACK #######
    drop <- 200
    n <- n+drop
    ####################
    
    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()
      }
    }
    
    # modifs : on retourne l'information d'arrivée, de départ
    # t2-t1
    res <- data.frame()
    ###### HACK #######
    #for ( i in 1:n ){
    for ( i in drop:n ){
    ###################
      #res=rbind(res,c(A[i], t1[i],t2[i],S[i], t2[i]-t1[i], t2[i]-t1[i]-S[i]))
      res=rbind(res, t2[i]-t1[i])
    }
    #colnames(res)<-c("inter", "arr","dep", "serv", "tps", "att")
    colnames(res)<-c("tps")
    return(res)
}    

Nous réalisons quelques tests préalables pour prendre en main le problème et pour vérifier que nos modifications font ce que nous souhaitons sur de petits exemples :

FileLIFO(n=5,lambda=1, typeservice="det",x=0, y=1,policy="npmtn" )
##        tps
## 1 1.605459
## 2 5.793085
## 3 4.663440
## 4 1.415016
## 5 2.970198
## 6 1.468726
FileLIFO(n=5,lambda=10,typeservice="uni",x=0, y=1,policy="npmtn" )
##         tps
## 1 3.0600540
## 2 2.4390589
## 3 1.4931001
## 4 0.9491297
## 5 0.5434993
## 6 0.6075883
FileLIFO(n=5,lambda=10,typeservice="gamma",x=1, y=1,policy="npmtn" )
##        tps
## 1 9.081501
## 2 8.392283
## 3 7.937064
## 4 6.888289
## 5 4.229046
## 6 1.379219

Question 1 : Nature des lois de service

Dans cette partie, nous nous intéressons à jouer avec la fonction Service et de nous rendre compte des comportements des différentes lois.

Service prend 4 paramètres

Nous cherchons aussi ici à déterminer les paramètres de chacune de ces lois pour obtenir une espérance de 1 et donc être capables par la suite d’évaluer l’impact des lois de service sur les performances du système.

det

Il s’agit de la loi la plus simple des quatre lois : Elle retourne simplement un vecteur de n 1. Elle a donc bien une espérance de 1.

uni

Il s’agit d’une loi uniforme entre les valeurs x et y.

Nous réalisons une fonction qui permet de tracer les résultats de la loi

draw_function <- function(data, log=F,title="values generated by the law"){

  # draw results
  g <- ggplot(data, aes(x=i, y=time))+geom_bar(data=data,stat="identity", position="identity")
  
  if ( log )
    g <- g + scale_y_log10()
  
              
  # draw red line
  g <- g + geom_abline(intercept=mean(data$time),slope=0, colour="red")
  
  g <- g + ggtitle(title)
  
  g
}
draw_uni <- function(n,min,max){
  data <- data.frame(time=Service(n,"uni",min,max))
  data$i <- 1:n
  draw_function(data)
}

N <- 100
draw_uni(n=N,min=0,max=2)

Nous remarquons que nous sommes très proches de ce que nous attendions.

Cherchons maintenant à avoir une valeur pour min et max telle que l’espérance soit de 1 :

Nous prenons min=0 et max=2 et avons donc une espérance de maxmin2=22=1

gamma

Dans le cas de la loi Gamma, x correspond au paramètre shape et y a scale.

Testons l’influence de ces paramètres en réalisant une fonction draw_gamma qui prend en paramètre

  • n : nombre d’échantillons
  • x : vector de shape
  • y : vector de scale
  • nbin : nombre de bins pour l’histogramme
  • et qui retourne le plot correspondant

et trace autant de graphes que de valeurs dans x ou y.

draw_gamma <- function(n,x,y,nbin){
  # generating data
  data <- data.frame(graph=paste("x:",x[1]," - y:",y[1]), time=Service(n,"gamma",x[1],y[1]))
  for ( i in 2:length(x) ){
    data <- rbind(data, data.frame(graph=paste("x:",x[i]," - y:",y[i]), time=Service(n,"gamma",x[i],y[i])))
  }

  # draw histograms
  width <- (max(data$time)-min(data$time)+0.05)/nbin
  g <- ggplot(data, aes(x=time,fill=..count..))+geom_histogram(binwidth=width )
  g <- g + facet_grid(graph~.)

  g
}

draw_gamma(n=10000,x=c(1,2,4,8),y=c(1,1,1,1),nbin=100)

En faisant varier le paramètre x, on remarque que le paramètre x vait varier la valeur des valeurs les plus courantes.

draw_gamma(n=10000,y=c(1,2,4,8),x=c(1,1,1,1),nbin=100)

En faisant varier le paramètre y, on remarque que celui-ci ne fait varier que la variance.

Essayons maintenant de trouver deux valeurs de x et y telle que l’espérance de la loi gamma vaille 1 :

Nous fixons le paramètre x à 2 puis à 4 et cherchons à trouver le y qui nous donne une loi d’espérance 1. Nous cherchons le paramètre y par dichotomie avec l’expérimentation suivante :

x = 2
y = 0.5
n = 10000000
mean(Service(n,"gamma",x,y))
## [1] 1.000089

Nous en déduisons donc les valeurs de notre premier couple de paramètres qui garantit une espérance de 1 : (x=2,y=0.5)

Affichons quelques valeurs:

x = 2
y = 0.5
n = 100
data <- data.frame(time=Service(n,"gamma",x,y))
data$i <- 1:n
draw_function(data)

Refaisons l’expérience pour x=4

x = 4
y = 0.25
n = 10000000
mean(Service(n,"gamma",x,y))
## [1] 0.9996478

Nous avons donc notre deuxième couple de paramètres (x=4,y=0.25)

Affichons quelques valeurs :

x = 4
y = 0.25
n = 100
data <- data.frame(time=Service(n,"gamma",x,y))
data$i <- 1:n
draw_function(data)

Nous pouvons maintenant nous intéresser à l’écart type de ces différentes lois :

n = 100000
sd(Service(n,"gamma",2,0.5))
## [1] 0.7088033
sd(Service(n,"gamma",4,0.25))
## [1] 0.4986812

exp

Cette loi prend seulement un paramètre x.

Étudions celle ci en dessinant quelques graphes :

draw_exp <- function(n,x,nbin){
  # generating data
  data <- data.frame(graph=paste("x:",x[1]), time=Service(n,"exp",x[1],0))
  for ( i in 2:length(x) ){
    data <- rbind(data, data.frame(graph=paste("x:",x[i]), time=Service(n,"exp",x[i],0)))
  }

  # draw histograms
  width <- (max(data$time)-min(data$time)+0.05)/nbin
  g <- ggplot(data, aes(x=time,fill=..count..))+geom_histogram(binwidth=width )
  g <- g + facet_grid(graph~.)

  g
}

draw_exp(n=5000,x=c(0.5,1,2,4),nbin=100)

On note que plus la valeur de x est élevée, plus les valeurs générées seront proches de 0.

Cherchons la valeur de x qui donne une espérance égale à 1 :

x = 1
y = 0
n = 10000000
mean(Service(n,"exp",x,y))
## [1] 0.9994748

Nous obtenons x=1.

Traçons quelques valeurs :

x = 1
n = 100
data <- data.frame(time=Service(n,"exp",x,0))
data$i <- 1:n
draw_function(data)

Cette loi a pour écart type :

x = 1
n = 10000000
sd(Service(n,"exp",x,0))
## [1] 0.9997963

Conclusions

Dans cette question, nous avons pu trouver les valeurs des paramètres x et y pour avoir des lois d’espérance 1. Cela nous permettra par la suite de comparer les résultats en fonction des lois de service.

Rappelons les écarts type obtenues pour les différentes lois

n <- 10000
sd(Service(n,"det",0,0))
## [1] 0
sd(Service(n,"uni",0,2))
## [1] 0.5788687
sd(Service(n,"exp",1,0))
## [1] 0.9885149
sd(Service(n,"gamma",2,0.5))
## [1] 0.7139836
sd(Service(n,"gamma",4,0.25))
## [1] 0.5085317

Question 2 : Étude détaillée de la file M/M/1 - LIFO

Nous définissons une surcouche de FileLIFO qui choisit le paramètre typeservice pour être une file LIFO M/M/1.

MM1 <- function(n, lambda, x, policy){ 
  res <- FileLIFO(n,lambda,"exp", x, 0, policy)
  res$lambda=lambda
  return(res[c("tps","lambda")])
}

Intéressons nous maintenant à tracer l’évolution du temps de service en fonction de λ et des politiques. Nous définissons donc une fonction MM1_all_policies qui pour un n, lambda et x fixé renvoie une data frame avec l’ensemble des mesures pour chaque politique.

MM1_all_policies <- function(n, lambda, x){
  res <- MM1(n,lambda,x,"npmtn")
  res$policy <- "npmtn"
  l <- c("pmtn","pmtn_restart","pmtn_reset")
  for ( i in l ){
    tmp <- MM1(n,lambda,x,i)
    tmp$policy<-i
    res<-rbind(res,tmp)
  }
  return(res)
}

De la même façon que précédemment, nous définissons une fonction gen_MM1 qui génère les résultats pour chaque lambda avec les paramètres n et x.

gen_MM1 <- function(n,x){
  # generating data
  data <- MM1_all_policies(n,0.2,x)
  l <- c(0.4, 0.6, 0.8)
  for ( i in l ){
    data <- rbind(data,MM1_all_policies(n,i,x))
  }
  data
}

On en profite pour fixer un nombre N pour notre nombre d’éléments ainsi que la précision de nos intervalles de confiance:

N <- 2500
confidence <- 0.975

Nous pouvons maintenant réaliser une fonction draw_MM1 qui affiche les courbes de temps de service pour un n et x fixé. Nous donnons aussi la possibilité d’avoir les ordonnées avec une échelle logarithmique dans le cas où nous avons un outsider.

draw_MM1 <- function(n,x, log=F, title="temps moyen passé dans le système dans une file M/M/1 LIFO"){
  
  data <- gen_MM1(n,x)
  
  data <- ddply(data, c("lambda","policy"), summarise,
            mean=mean(tps),
            sd=sd(tps)
          )
  data$diff <- qnorm(confidence)*data$sd/sqrt(n)
  
  # draw graphs
  g <-ggplot(data, aes(x=lambda, y=mean, color=policy))+geom_line()
  if ( log )
    g <- g + scale_y_log10()
  
  # draw confidence intervals
  g <- g + geom_errorbar(width=0.01, aes(x=lambda, y=mean, ymin=mean-diff, ymax=mean+diff))
  
  g <- g + ggtitle(title)
  
  g
}
draw_MM1(N,1,log=T)