Simulation

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
set.seed(100)

mu1 = 1
mu2 = 2
veille1 = 0.5
veille2=0.25
non_veille1 = 1
non_veille2 = 4

simul = function(N=100, arrival_rate, processing_rate, debug=FALSE) {
  Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
  Service = rexp(n=N, rate =processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)

  Remaining = rep(N, x=NA);
  Completion = rep(N, x=NA)
  t = 0;
  
  
  CurrentTask = NA;
  NextArrival = 1;
  while (TRUE) {
    if(debug) print(t)
    if(debug) print(Arrival);
    if(debug) print(Service);
    if(debug) print(Remaining);
    dtA = NA;
    dtC = NA;
    
    
    if(length(Arrival[Arrival>t])>0) {
      dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
    }
    if(!is.na(CurrentTask)) {
      dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
    }
    if(is.na(dtA) & is.na(dtC)) {
      break;
    } 
    dt = min(dtA,dtC,na.rm=T)
    
    # Mettre à jour comme il faut:
    #   la date
    t = t + dt;
    #   les arrivées
    if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
      Remaining[NextArrival] = Service[NextArrival];
      NextArrival = NextArrival + 1;
    }
    
    #   le remaining 
    if(!is.na(CurrentTask)) {
      Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
      if(Remaining[CurrentTask] <= 0) {
        Completion[CurrentTask] = t;
        Remaining[CurrentTask] = NA;
      }
      CurrentTask = NA;
    }
    #   et currentTask (politique d'ordonnancement: FIFO)
    WaitingList=(1:N)[!is.na(Remaining)];
    if(length(WaitingList)>0) {
      CurrentTask = head(WaitingList,n=1);
    }
    
  }
  
  if(debug) print(Completion)
  return(data.frame(
    Arrival_rate=arrival_rate,
    Processing_rate=processing_rate,
    Arrival = Arrival,
    Service = Service,
    Completion = Completion
  ))
  
}

Question 1 : Comparaison de la performance de S1 et de S2

Espérance du temps de réponse des deux serveurs

temps_de_reponse_1  = function(min,max,step){
  df = data.frame()
  vec=seq(from=min, to=max, by=step)
  for (r in  vec) {
    df = rbind(df,simul(N=1000, arrival_rate=r, processing_rate=mu1, debug=FALSE));
  }
  
  return(df) 
  
}
#####################################################################################
temps_de_reponse_2 = function(min,max,step){
  df4 = data.frame()
  vec=seq(from=min, to=max, by=step)
  for (r in  vec) {
    df4 = rbind(df4,simul(N=1000, arrival_rate=r, processing_rate=mu2, debug=FALSE));
  }
  
  return(df4) 
  
}
######################################################################################""
#set.seed(100);
# on calcule le temps de réponse de s1 et s2 
min = 0.1
max = 0.9
step = 0.1

df4 = temps_de_reponse_2(min,max,step)
df = temps_de_reponse_1(min,max,step)


df2 = df %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response = mean(Completion - Arrival), 
            Response_se = sd(Completion - Arrival)/sqrt(n()))

df5 = df4 %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response1 = mean(Completion - Arrival), 
            Response_se1 = sd(Completion - Arrival)/sqrt(n()))

df_complet=cbind(df2,df5)
p = ggplot(df_complet, aes(x=Arrival_rate, y=Response)) + geom_line(color=2)  +
  geom_errorbar(aes(ymin = Response - 2*Response_se,
                    ymax = Response + 2*Response_se), color=2)


p = p + geom_line(data=df_complet, aes(x=Arrival_rate, y=Response1))  +
 geom_errorbar(aes(ymin = Response1 - 2*Response_se1,
                  ymax = Response1 + 2*Response_se1))

p = p + ggtitle("       les temps de réponse de s1 et s2 en fonction de lambda")
p

# on augmente les valeurs de lambda pour voir un peu plus claire les courbes 
min = 0.1
max = 1.9
step = 0.1

df4 = temps_de_reponse_2(min,max,step)
df = temps_de_reponse_1(min,max,step)

df2 = df %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response = mean(Completion - Arrival), 
            Response_se = sd(Completion - Arrival)/sqrt(n()))

df5 = df4 %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response1 = mean(Completion - Arrival), 
            Response_se1 = sd(Completion - Arrival)/sqrt(n()))

df_complet=cbind(df2,df5)
p = ggplot(df_complet, aes(x=Arrival_rate, y=Response)) + geom_line(color=2) +  ylim(0,50)+
  geom_errorbar(aes(ymin = Response - 2*Response_se,
                    ymax = Response + 2*Response_se), color=2)


p = p + geom_line(data=df_complet, aes(x=Arrival_rate, y=Response1))  +
 geom_errorbar(aes(ymin = Response1 - 2*Response_se1,
                  ymax = Response1 + 2*Response_se1))

p = p + ggtitle("       les temps de réponse de s1 et s2 en fonction de lambda ")

p
## Warning: Removed 9 rows containing missing values (geom_path).
## Warning: Removed 9 rows containing missing values (geom_errorbar).

Espérance de la consommation énérgétique des deux serveurs

simul_1serveur = function(N=100, arrival_rate=.02, processing_rate, debug=FALSE, veille , non_veille) {
  Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
  Service = rexp(n=N, rate =processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)

  Remaining = rep(N, x=NA);
  Completion = rep(N, x=NA)
  t = 0;
  consom = 0 
  
  CurrentTask = NA;
  NextArrival = 1;
  while (TRUE) {
    if(debug) print(t)
    if(debug) print(Arrival);
    if(debug) print(Service);
    if(debug) print(Remaining);
    dtA = NA;
    dtC = NA;
    
    
    if(length(Arrival[Arrival>t])>0) {
      dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
    }
    if(!is.na(CurrentTask)) {
      dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
    }
    if(is.na(dtA) & is.na(dtC)) {
      break;
    } 
    dt = min(dtA,dtC,na.rm=T)
    
    # Mettre à jour comme il faut:
    #   la date
    t = t + dt;
    #   les arrivées
    if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
      Remaining[NextArrival] = Service[NextArrival];
      NextArrival = NextArrival + 1;
    }
    
    #   le remaining 
    if(!is.na(CurrentTask)) {
      consom = consom + non_veille * dt
      Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
      if(Remaining[CurrentTask] <= 0) {
        Completion[CurrentTask] = t;
        Remaining[CurrentTask] = NA;
      }
      CurrentTask = NA;
    }
    
    else{
      consom = consom + veille * dt
      
    }
    #   et currentTask (politique d'ordonnancement: FIFO)
    WaitingList=(1:N)[!is.na(Remaining)];
    if(length(WaitingList)>0) {
      CurrentTask = head(WaitingList,n=1);
    }
    
  }
  
  if(debug) print(Completion)
  return(data.frame(
    Arrival_rate=arrival_rate,
    Processing_rate=processing_rate,
    Arrival = Arrival,
    consom = consom/t , 
    Completion = Completion
  ))
  
}


df = data.frame()

for (r in  seq(.1,.2,.1)) {
  df = rbind(df, data.frame(simul_1serveur(N = 1000, arrival_rate = r, 1,FALSE,0.5,1),s="1"))
}

for (r in  seq(.1,.2,.1)) {
  df = rbind(df, data.frame(simul_1serveur(N = 1000, arrival_rate = r,processing_rate = 2,FALSE,0.25,4), s ="2"))
}

df2 = df %>% group_by(Arrival_rate, s) %>%
    summarize(consom = mean(consom))


ggplot(df2, aes(x=Arrival_rate, y=consom, group = s, color = s)) + geom_line() 

En zoomant sur le point d’intersection on trouve une valeur approché de \(\lambda = 0.18\) #Rentabilité énérgétique Pour λ>0.19 il est donc préférable en terme de consommation de choisir le serveur1 plutôt que le serveur2. Cela semble être cohérent parce que le serveur2 est concu pour etre performant mais energivore quand il s’agit de taux d’arrivée elevé, S2 consomme 4 alors que S1 en consomme 1. Pour des taux d’arrivée faible, le S2 fini de traiter ses tâches rapidement, donc il est souvent en veille, S2 consomme 2 fois moins.

Question 2 : Etude du serveur S3

On a maintenant deux serveurs différents (chacun leur \(\lambda\) et leur \(\mu\)) qui travaillent en parallèle.

# simulation de 2 serveurs qui travaillent en parralléle 
# argument de simul + mu2 + temps de veille et de non_veille de 2éme serveur 
# retour , dataframe + consommation moyenne 
simul2 = function(N=100, arrival_rate=.02, processing_rate=1 ,processing_rate2=2 ,veille1 , non_veille1 , veille2 , non_veille2 ,debug=FALSE) {
  Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
  Service1 = rexp(n=N, rate =processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
  Service2 = rexp(n=N, rate =processing_rate2); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
   
  Remaining = rep(N, x=NA);
  Remaining2 = rep(N, x=NA)
  Completion = rep(N, x=NA)
  t = 0;
  consom = 0
  CurrentTask = NA;
  CurrentTask2 = NA;
  NextArrival = 1;
  #consom[i] = veille1 * Arrival[i]
  while (TRUE) {
    if(debug) print(t)
    if(debug) print(Arrival);
    if(debug) print(Service);
    if(debug) print(Remaining);
    dtA = NA;
    dtC = NA;
    dtC2 = NA
    if(length(Arrival[Arrival>t])>0) {
      dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
    }
    if(!is.na(CurrentTask)) {
      dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
    }
    
    if(!is.na(CurrentTask2)) {
      dtC2 = Remaining2[CurrentTask2]; # temps jusqu'à la prochaine terminaison
    }
    if(is.na(dtA) & is.na(dtC)  & is.na(dtC2)) {
      break;
    } 
    dt = min(dtA,dtC,dtC2,na.rm=T)

    # Mettre à jour comme il faut:
    #   la date
    t = t + dt;
    
    
    #   les arrivées
    if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
    
      Remaining[NextArrival] = Service1[NextArrival]
      Remaining2[NextArrival] = Service2[NextArrival]
      NextArrival = NextArrival + 1;
    }
    
   
    
    #   le remaining 
    if(!is.na(CurrentTask)) {
      Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
      if(Remaining[CurrentTask] <= 0) {
        Completion[CurrentTask] = t;
        Remaining[CurrentTask] = NA;
        
      }
      consom = consom + non_veille1 * dt 
      
      CurrentTask = NA;
    }
    else {
      consom = consom+   veille1 * dt
    }
    
    if(!is.na(CurrentTask2)) {
      
      Remaining2[CurrentTask2] = Remaining2[CurrentTask2] - dt ;
      if(Remaining2[CurrentTask2] <= 0) {
        Completion[CurrentTask2] = t;
        Remaining2[CurrentTask2] = NA;
        
      }
      consom =  consom +non_veille2 * dt
      CurrentTask2 = NA;
    }
    else{
      consom =  consom +veille2 * dt
        
    }
    #   et currentTask (politique d'ordonnancement: FIFO)
    WaitingList=(1:N)[!is.na(Remaining)];
    if(length(WaitingList)>0) {
      
       CurrentTask = head(WaitingList,n=1);
       Remaining2[CurrentTask] = NA
       
    }
    
    WaitingList2=(1:N)[!is.na(Remaining2)];
    if(length(WaitingList2)>0) {
      CurrentTask2 = head(WaitingList2,n=1);
      Remaining[CurrentTask2] = NA
    }
    
  }
  
  if(debug) print(Completion)
  return(data.frame(
    Arrival_rate=arrival_rate,
    Processing_rate=processing_rate,
    Arrival = Arrival,
    Completion = Completion,
    consom = consom/t 
 
  ))
  
}




# calculer temps de réponse de serveur 3 
temps_de_reponse_3 = function(min,max,step){
  df = data.frame()
  vec=seq(from=min, to=max, by=step)
  for (r in  vec) {
    df = rbind(df,simul2(N=100, arrival_rate=r,mu1,mu2,veille1,non_veille1,veille2,non_veille2,debug=FALSE));
  }
  
  return(df) 
  
}
min = 0.1
max = 2
step = 0.1


# calculer temps de réponse de s2 et s3 et les dessiner 
df = data.frame()
df4= data.frame()
df = temps_de_reponse_3(min,max,step)
df4 = temps_de_reponse_1(min,max,step)
          
df2 = df %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response = mean(Completion - Arrival), 
            Response_se = sd(Completion - Arrival)/sqrt(n()))

df5 = df4 %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response1 = mean(Completion - Arrival), 
            Response_se1 = sd(Completion - Arrival)/sqrt(n()))



df_complet = data_frame()
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
df_complet=cbind(df2,df5)
p = ggplot(df_complet, aes(x=Arrival_rate, y=Response)) + geom_line()  +
  geom_errorbar(aes(ymin = Response - 2*Response_se,
                    ymax = Response + 2*Response_se), color=2)


p = p + geom_line(data=df_complet, aes(x=Arrival_rate, y=Response1))  +
  geom_errorbar(aes(ymin = Response1 - 2*Response_se1,
                    ymax = Response1 + 2*Response_se1)) 

p = p + ggtitle("       la temps de reponse de s2 et s3 en fonction de lambda")

p

####consomation s2 et s3
# fonction comme la fonction classique simul mais on calcule la consommation avec 
set.seed(100)
df =data_frame()
for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul_1serveur(N = 1000, arrival_rate = r, 1,FALSE,0.5,1),s="1"))
}

for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul_1serveur(N = 1000, arrival_rate = r,processing_rate = 2,FALSE,0.25,4), s ="2"))
}


for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul2(N = 1000,arrival_rate = r, 1 , 2 ,debug =FALSE, 0.5 , 1 , 0.25 ,4), s ="3"))
}
df2 = df %>% group_by(Arrival_rate, s) %>%
    summarize(consom = mean(consom))


p = ggplot(df2, aes(x=Arrival_rate, y=consom, group = s, color = s)) + geom_line() 
p

On attend du S3 qu’il soit un compromis entre le temps de réponse et la consommation énergetique. Intuitivement, le S3 est équipé aussi bien pour des λ petits que des λ grands. D’après ce graphe, pour λ>0.60 il est légèrement préférable en terme de consommation de choisir le serveur3 plutôt que le serveur2. Cela semble être cohérent parce que le serveur3 est concu pour etre performant et moins energivore (que S2) quand il s’agit de taux d’arrivée elevé car les tâches seront traitées majoritairement par sa composante de type S2 et de temps en temps par son autre composante moins performant mais moins energivore. Cette dernière agit en effet comme un adoucisseur de la consommation. Tandis que S2, pour des λ grands, devra travailler en continu, sans rien pour l’aider dans le traitement des tâches, S2 va donc être plus energivore.

Question 3 : Influence de l’hypothèse Markovienne

On suppose donc que le temps de service est déterministe, c’est a dire que toutes arrivée sont traitées pendant un temps de service égal. Nous voulons alors voir l’influence de ce caractère deterministe sur notre temps de réponse et notre consommation.

simul_markov_1serveur = function(N=100, arrival_rate, processing_rate, debug=FALSE , veille , non_veille) {
  Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
  #Service = rexp(n=N, rate =1/processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
  Service = rep(1/processing_rate,N);
  
  Remaining = rep(N, x=NA);
  Completion = rep(N, x=NA)
  t = 0;
  consom = 0 
  
  CurrentTask = NA;
  NextArrival = 1;
  while (TRUE) {
    if(debug) print(t)
    if(debug) print(Arrival);
    if(debug) print(Service);
    if(debug) print(Remaining);
    dtA = NA;
    dtC = NA;
    
    
    if(length(Arrival[Arrival>t])>0) {
      dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
    }
    if(!is.na(CurrentTask)) {
      dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
    }
    if(is.na(dtA) & is.na(dtC)) {
      break;
    } 
    dt = min(dtA,dtC,na.rm=T)
    
  
    t = t + dt;
    
    if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
      Remaining[NextArrival] = Service[NextArrival];
      NextArrival = NextArrival + 1;
    }
    
    if(!is.na(CurrentTask)) {
      Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
      if(Remaining[CurrentTask] <= 0) {
        Completion[CurrentTask] = t;
        Remaining[CurrentTask] = NA;
      }
      consom = consom + non_veille * dt
      CurrentTask = NA;
    }
    
    else{
      consom = consom + veille * dt
      
    }
    WaitingList=(1:N)[!is.na(Remaining)];
    if(length(WaitingList)>0) {
      CurrentTask = head(WaitingList,n=1);
    }
    
  }
  
  if(debug) print(Completion)
  return(data.frame(
    Arrival_rate=arrival_rate,
    Processing_rate=processing_rate,
    Arrival = Arrival,
    consom = consom/t , 
    Completion = Completion
  ))
  
}


temps_de_reponse_markov_1 = function(min,max,step){
  df = data.frame()
  vec=seq(from=min, to=max, by=step)
  for (r in  vec) {
    df = rbind(df,simul_markov_1serveur(N=1000, arrival_rate=r, processing_rate=mu1, debug=FALSE, veille1 , non_veille1));
  }
  
  return(df) 
  
}

temps_de_reponse_markov_2 = function(min,max,step){
  df4 = data.frame()
  vec=seq(from=min, to=max, by=step)
  for (r in  vec) {
    df4 = rbind(df4,simul_markov_1serveur(N=1000, arrival_rate=r, processing_rate=mu2, debug=FALSE , veille2,non_veille2));
  }
  
  return(df4) 
  
}

set.seed(100);
df = data.frame()
df4 = data.frame()
min = 0.1
max = 0.9
step = 0.1

df4 = temps_de_reponse_markov_2(min,max,step)
df = temps_de_reponse_markov_1(min,max,step)


df2 = df %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response = mean(Completion - Arrival), 
            Response_se = sd(Completion - Arrival)/sqrt(n()))

df5 = df4 %>% group_by(Arrival_rate, Processing_rate) %>%
  summarize(Response1 = mean(Completion - Arrival), 
            Response_se1 = sd(Completion - Arrival)/sqrt(n()))


df_complet=cbind(df2,df5)
p = ggplot(df_complet, aes(x=Arrival_rate, y=Response)) + geom_line(color=2)   +
  geom_errorbar(aes(ymin = Response - 2*Response_se,
                    ymax = Response + 2*Response_se), color=2)


p = p + geom_line(data=df_complet, aes(x=Arrival_rate, y=Response1))  +
 geom_errorbar(aes(ymin = Response1 - 2*Response_se1,
                  ymax = Response1 + 2*Response_se1))

p = p + ggtitle("       le temps de réponse de s1 et s2 en fonction de lambda 
                avec temps de service déterministe")


p

Pour la consommation :

df = data.frame()

for (r in  seq(.1,.2,.1)) {
  df = rbind(df, data.frame(simul_markov_1serveur(N = 1000, arrival_rate = r, 1,FALSE,0.5,1),s="1"))
}

for (r in  seq(.1,.2,.1)) {
  df = rbind(df, data.frame(simul_markov_1serveur(N = 1000, arrival_rate = r,processing_rate = 2,FALSE,0.25,4), s ="2"))
}

df2 = df %>% group_by(Arrival_rate, s) %>%
    summarize(consom = mean(consom))


ggplot(df2, aes(x=Arrival_rate, y=consom, group = s, color = s)) + geom_line()  

Les courbes de consommation ne change pas pour des serveurs travaillant seul

question 3.2

simul2_markov = function(N=100, arrival_rate=.02, processing_rate=1 ,processing_rate2=2 ,veille1 , non_veille1 , veille2 , non_veille2 ,debug=FALSE) {
  Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
  #Service1 = rexp(n=N, rate =1/processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
  #Service2 = rexp(n=N, rate =1/processing_rate2); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
  Service1 = rep(1/processing_rate,N);
  Service2 = rep(1/processing_rate2,N);
 
  
  Remaining = rep(N, x=NA);
  Remaining2 = rep(N, x=NA)
  Completion = rep(N, x=NA)
  t = 0;
  consom = 0
  CurrentTask = NA;
  CurrentTask2 = NA;
  NextArrival = 1;
  while (TRUE) {
    if(debug) print(t)
    if(debug) print(Arrival);
    if(debug) print(Service);
    if(debug) print(Remaining);
    dtA = NA;
    dtC = NA;
    dtC2 = NA
    if(length(Arrival[Arrival>t])>0) {
      dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
    }
    if(!is.na(CurrentTask)) {
      dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
    }
    
    if(!is.na(CurrentTask2)) {
      dtC2 = Remaining2[CurrentTask2]; # temps jusqu'à la prochaine terminaison
    }
    if(is.na(dtA) & is.na(dtC)  & is.na(dtC2)) {
      break;
    } 
    dt = min(dtA,dtC,dtC2,na.rm=T)

    # Mettre à jour comme il faut:
    #   la date
    t = t + dt;
    
    
    #   les arrivées
    if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
    
      Remaining[NextArrival] = Service1[NextArrival]
      Remaining2[NextArrival] = Service2[NextArrival]
      NextArrival = NextArrival + 1;
    }
    
   
    
    #   le remaining 
    if(!is.na(CurrentTask)) {
      Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
      if(Remaining[CurrentTask] <= 0) {
        Completion[CurrentTask] = t;
        Remaining[CurrentTask] = NA;
        
      }
      consom = consom + non_veille1 * dt 
      
      CurrentTask = NA;
    }
    else {
      consom = consom+   veille1 * dt
    }
    
    if(!is.na(CurrentTask2)) {
      
      Remaining2[CurrentTask2] = Remaining2[CurrentTask2] - dt ;
      if(Remaining2[CurrentTask2] <= 0) {
        Completion[CurrentTask2] = t;
        Remaining2[CurrentTask2] = NA;
        
      }
      consom =  consom +non_veille2 * dt
      CurrentTask2 = NA;
    }
    else{
      consom =  consom +veille2 * dt
        
    }
    WaitingList=(1:N)[!is.na(Remaining)];
    if(length(WaitingList)>0) {
      
       CurrentTask = head(WaitingList,n=1);
       Remaining2[CurrentTask] = NA
       
    }
    
    WaitingList2=(1:N)[!is.na(Remaining2)];
    if(length(WaitingList2)>0) {
      CurrentTask2 = head(WaitingList2,n=1);
      Remaining[CurrentTask2] = NA
    }
    
  }
  
  if(debug) print(Completion)
  return(data.frame(
    Arrival_rate=arrival_rate,
    Processing_rate=processing_rate,
    Arrival = Arrival,
    Completion = Completion,
    consom = consom/t 
    
  ))
  
}

df = data.frame()


for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul_markov_1serveur(N = 1000, arrival_rate = r,processing_rate = 2,FALSE,0.25,4), s ="2"))
}
for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul2_markov(N = 1000, arrival_rate = r,1,2,0.5,1,0.25,2,FALSE), s ="3"))
}


df2 = df %>% group_by(Arrival_rate, s) %>%
    summarize(Response = mean(Completion - Arrival), 
              Response_se = sd(Completion - Arrival)/sqrt(n()))

ggplot(df2, aes(x=Arrival_rate, y=Response, group = s, color = s)) + geom_line()   

Pour 2 serveurs travaillant en parallèle. On observe quelque chose de bizarre, la courbe est en effet decroissante et ne se croise pas avec la courbe de s2. C’est là ou on remarque l’importance de l’hypothese markovienne.

energie markov s2 vs s3

set.seed(100)
df =data_frame()
for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul_1serveur(N = 1000, arrival_rate = r, 1,FALSE,0.5,1),s="1"))
}

for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul_1serveur(N = 1000, arrival_rate = r,processing_rate = 2,FALSE,0.25,4), s ="2"))
}


for (r in  seq(.1,.9,.1)) {
  df = rbind(df, data.frame(simul2_markov(N = 1000,arrival_rate = r, 1 , 2 ,debug =FALSE, 0.5 , 1 , 0.25 ,4), s ="3"))
}
df2 = df %>% group_by(Arrival_rate, s) %>%
    summarize(consom = mean(consom))


p = ggplot(df2, aes(x=Arrival_rate, y=consom, group = s, color = s)) + geom_line()   
p

On constate que les graphiques obtenus sont très similaires à ceux obtenus avec des temps de service exponentiels . c’est normal vu que l’espérence de loi exponontielle de paramétre lambda = 1/lambda .