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 :

• inter : temps d’interarrivée
• arr : date d’arrivée
• dep : date de départ
• serv : temps de service
• tps : temps passé dans le système ($arr-dep$)
• att : temps passé à attendre dans le système ($tps-serv$, dans le meilleur des cas, vaut 0)

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

• n : qui caractérise la taille du vecteur réponse
• typeservice : qui permet de choisir la loi a utiliser
• x et y : la signification dépendent de la loi utilisée

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.

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)

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)

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

x = 2
y = 0.5
n = 10000000
mean(Service(n,"gamma",x,y))

## [1] 1.000089

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

x = 4
y = 0.25
n = 10000000
mean(Service(n,"gamma",x,y))

## [1] 0.9996478

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

n = 100000
sd(Service(n,"gamma",2,0.5))

## [1] 0.7088033

sd(Service(n,"gamma",4,0.25))

## [1] 0.4986812

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)

x = 1
y = 0
n = 10000000
mean(Service(n,"exp",x,y))

## [1] 0.9994748

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

x = 1
n = 10000000
sd(Service(n,"exp",x,0))

## [1] 0.9997963

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 $\lambda$ 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)

Pour cette courbe, nous utilisons l’échelle logarithmique car restart est beaucoup plus élevé que les autres courbes. Nous remarquons que les autres politiques donnent des résultats similaires

Sur l’expérience que nous avons réalisée, nous remarquons que la politique restart est bien moins bonne que les autres. Nous pensons que cela vient du fait que comme le processus est interrompu et redémarré quand un autre arrive, il pourra rester longtemps si son temps de service est important car il aura de fortes chances de se faire doubler régulièrement avant de pouvoir sortir.

Nous pensons que la politique reset s’en sort mieux que le restart car un processus qui s’est fait doubler prend un nouveau temps de service, qui peut être petit et donc, il pourra sortir plus vite de la file.

Question 3 : Étude de la file M/GI/1 - LIFO

Dans cette partie, nous nous intéressons à refaire les expérimentations de la question 2, mais avec d’autres loi de temps de service.

Dans un premier temps, nous définissons des fonctions qui renvoient les résultats de simulation avec les différentes lois de temps de service.

Simulation avec un temps de service constant :

cst <- function(n, lambda,x,y, policy){
  res <- FileLIFO(n,lambda,"det",0,0,policy)
  return(res["tps"])
}

Simulation avec un temps de service de loi uniforme :

uni <- function(n,lambda,x,y,policy){
  res <- FileLIFO(n,lambda,"uni",x,y,policy)
  return(res["tps"])
}

Simulation avec une loi gamma :

gamma1 <- function(n,lambda,x,y,policy){
  res <- FileLIFO(n,lambda,"gamma",2,0.5,policy)
  return(res["tps"])
}

gamma2 <- function(n,lambda,x,y,policy){
  res <- FileLIFO(n,lambda,"gamma",4,.25,policy)
  return(res["tps"])
}

Nous pouvons maintenant réaliser une fonction qui génère les données du graphe en fonction de la fonction donnée en paramètre, de x et de y.

Nous définissons deux fonctions :

• generate_data_lambda qui génère tous les résultats pour la fonction f et pour les paramètres lambda, x, y.

• generate_data qui reprend generate_data_lambda mais en générant pour les lambdas 0.2, 0.4, 0.6 et 0.8.

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

generate_data <- function(f, n, x, y){
  res <- generate_data_lambda(f,n,0.2,x,y)
  res<-rbind(res,generate_data_lambda(f,n,0.4,x,y))
  res<-rbind(res,generate_data_lambda(f,n,0.6,x,y))
  res<-rbind(res,generate_data_lambda(f,n,0.8,x,y))
  return(res)
}

Et enfin, une fonction affichant les données ainsi générées:

draw_data <- function(data,n,log=F, title="temps passé dans le système"){
  
  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
  
}

data <- generate_data(cst, N, 0, 0)
draw_data(data,N, log=T, title="temps passé dans le système pour une loi de service déterministe")

data <- generate_data(uni, N, 0, 2)
draw_data(data,N,log=T, title="temps passé dans le système pour une loi de service uniforme")

data <- generate_data(gamma1, N, 0, 0)
draw_data(data,N,log=T, title="temps passé dans le système pour une loi de service gamma1")

data <- generate_data(gamma2, N, 0, 0)
draw_data(data,N,log=T, title="temps passé dans le système pour une loi de service gamma2")

Conclusion

Forts de ces expériences, nous remarquons que la politique pmtn_restart est moins bonne que les autres de manière générale et que npmtn est la meilleure suivie de pmtn.

• npmtn ressemble à une FIFO car le processus qui rentre s’exécute jusqu’au bout et ce sera le dernier entré qui aura accès en premier. Cependant, comme les processus ont des temps de services indépendants les uns des autres, le système se comporte comme une FIFO.

• pmtn est la politique préemptive simple, c’est à dire qui ne redémarre pas l’exécution des processus quand ils se font doubler.

Sur ces expériences, nous remarquons que la loi de service n’influe que peu sur les politiques à l’exception de pmtn_reset. On remarque que plus la loi de service à une variance importante, mieux pmtn_reset à une performance qui se rapproche des autres (cas de la loi exponentielle), plus la variance est faible, plus elle se rapproche de pmtn_restart (cas de la loi déterministe).

Cela peut s’expliquer par le fait que quand la variance tend vers 0, le fait de tirer une nouvelle durée de service ne change pas le temps de service. la politique reset est donc similaire à restart dans ce cas.

Si la variance est importante, quand une longue durée est tirée, l’espérance étant de 1, il y a de fortes chances qu’une plus petite durée soit tirée la prochaine fois ce qui permettra à la tâche de se terminer rapidement. Dans le cas du reset, la tâche risque d’être bloquée jusqu’à ce qu’il n’y ait plus de tâches dans le système.

Ce DM nous permet donc de dire que si nous avons une politique restart ou si nous avons une politique reset avec une variance faible, le système risque d’avoir de mauvaises performances tandis que les politiques npmtn et pmtn ont de bonnes performances quelle que soit la loi de service.