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 FIFO. Le processus d'arrivée est un processus de Poisson de taux \lambda (débit), les clients ont un temps de service de moyenne 1 pris comme unité de temps de référence.
Pour simuler une file FIFO, il suffit d'utiliser la méthode équationnelle à partir de la relation de récurrence sur les temps de réponse (équation de /Lindley/)
library(plyr)
library(ggplot2)
Service <- function(n = 1, typeservice, x, y) {
# genere un temps de service
switch(typeservice, det = seq(from = 1, to = 1, length.out = n), uni = runif(n,
x, y), gamma = rgamma(n, shape = x, scale = y), exp = rexp(n, x))
}
Response <- function(n, lambda, typeservice, x, y) {
# genere une trajectoire de temps de reponse pour une file M/GI/1 FIFO
# parametres : n : taille de la trajectoire lambda : debit d'arrivee
# typeservice : type de service (det uni gamma exp) x ,y : parametres de la
# loi de service valeur de retour : vecteur de n valeurs successives du
# temps de réponse dans une M/GI/1 FIFO initialement vide
# generation des arrivees et des temps de service
A <- rexp(n, lambda)
t1 <- cumsum(A)
S <- Service(n, typeservice, x, y)
# initialisation des structures (pour éviter une allocation dynamique et
# gagner du temps)
t2 <- seq(0, 0, length.out = n)
# Par curiosité comptons également le nombre de fois où le système se vide
sys_vide = 1
# initialisation du premier terme
t2[1] <- t1[1] + S[1]
# calcul de la trajectoire equation de recurrence de Lindley
for (i in 2:n) {
t2[i] <- S[i] + max(t1[i], t2[i - 1])
if (t1[i] > t2[i - 1]) {
sys_vide = sys_vide + 1
}
}
# rajoutons à la sortie, par curiosité, l'écart-type du temps de service...
list(response = t2 - t1, sys_vide = sys_vide, sd_S = sd(S))
}
create <- function(n, lambda_min, lambda_max, step, typeservice, x, y) {
d <- data.frame()
for (lambda in seq(lambda_min, lambda_max, step)) {
traj = Response(n, lambda, typeservice, x, y)
d = rbind(d, data.frame(n = n, lambda = lambda, response = mean(traj$response),
sd = sd(traj$response), sys_vide = traj$sys_vide, sd_S = traj$sd_S))
}
d$input_type = typeservice
d$input_p1 = x
d$input_p2 = y
d$label = as.factor(paste(typeservice, "(", x, ",", y, ")", sep = ""))
d
}
Pour chaque loi et pour des valeurs de \( \lambda \) de 0.01 à 0.95 avec un pas de 0.1 on estime le temps de réponse moyen des clients sur des échantillons de taille \( N=10000 \).
experiment <- function(sample_size = 10000) {
resp_exp <- create(sample_size, 0.01, 0.95, 0.05, "exp", 1, 0)
resp_det <- create(sample_size, 0.01, 0.95, 0.05, "det", 0, 0)
resp_u02 <- create(sample_size, 0.01, 0.95, 0.05, "uni", 0, 2)
resp_u0.5_1.5 <- create(sample_size, 0.01, 0.95, 0.05, "uni", 0.5, 1.5)
resp_gamma4 <- create(sample_size, 0.01, 0.95, 0.05, "gamma", 4, 1/4)
resp_gammap2 <- create(sample_size, 0.01, 0.95, 0.05, "gamma", 0.2, 1/0.2)
rbind(resp_exp, resp_det, resp_u02, resp_u0.5_1.5, resp_gamma4, resp_gammap2)
}
df <- experiment()
head(df)
## n lambda response sd sys_vide sd_S input_type input_p1 input_p2
## 1 10000 0.01 1.017 1.034 9901 1.0226 exp 1 0
## 2 10000 0.06 1.053 1.037 9357 0.9918 exp 1 0
## 3 10000 0.11 1.115 1.107 8916 0.9928 exp 1 0
## 4 10000 0.16 1.180 1.163 8429 1.0050 exp 1 0
## 5 10000 0.21 1.299 1.294 7877 1.0350 exp 1 0
## 6 10000 0.26 1.347 1.377 7444 1.0143 exp 1 0
## label
## 1 exp(1,0)
## 2 exp(1,0)
## 3 exp(1,0)
## 4 exp(1,0)
## 5 exp(1,0)
## 6 exp(1,0)
Regardons maintenant comment évolue le temps de réponse en fonction du taux d'arrivée et de la loi du temps de service:
ggplot(df, aes(x = lambda, y = response, color = label, shape = label)) + geom_line() +
geom_point() + geom_errorbar(width = 0.02, aes(x = lambda, y = response,
ymin = response - 2 * sd/sqrt(n), ymax = response + 2 * sd/sqrt(n))) + geom_vline(xintercept = 1) +
geom_hline(yintercept = 1) + theme_bw()
Mmmh, bon, on a des courbes qui “remontent” et qui “descendent”“ malgré des intervalles de confiance relativement petits. En fait, les temps de réponses sont correllés les uns aux autres et donc, on ne peut pas calculer les intervalles de confiance comme ça. Par contre, quand le système se vide, c'est bon, on a bien une remise à 0 et les temps de réponse redeviennent indépendant de ce qui s'est passé avant. Pareil, il faudrait prendre la variance de la moyenne temps de réponse entre deux resets mais là, j'ai un peu la flemme donc restons avec le code tel qu'il est et utilisons sys_vide plutôt que n pour calculer notre intervalle de confiance.
D'autre part, quand le système devient très très instable, si le temps de réponse augmente comme on pouvait s'y attendre, ça n'a pas grand sens de considérer le cas où le système a des temps de réponse super longs. En effet, dans un vrai système, quand le temps de réponse devient trop important, d'autres phénomènes apparaissent (les clients quittent le système, la file d'attente est tellement longue que le système se met à swapper, …). On va donc se concentrer sur le cas où le temps de réponse reste "raisonnable”.
ggplot(df, aes(x = lambda, y = response, color = label)) + geom_line(aes(linetype = input_type)) +
geom_point() + geom_errorbar(width = 0.02, aes(x = lambda, y = response,
ymin = response - 2 * sd/sqrt(sys_vide), ymax = response + 2 * sd/sqrt(sys_vide))) +
geom_vline(xintercept = 1) + geom_hline(yintercept = 1) + theme_bw() + scale_color_brewer(palette = "Dark2") +
ylim(0, 6) + 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")
## Warning: Removed 10 rows containing missing values (geom_path).
## Warning: Removed 10 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 4 rows containing missing values (geom_path).
Bon, au final, on a toujours des intervalles de confiance un peu large. Donc, recommençons avec une trajectoire plus importante pour lisser les choses.
df <- experiment(sample_size = 1e+05)
head(df)
## n lambda response sd sys_vide sd_S input_type input_p1 input_p2
## 1 1e+05 0.01 1.007 1.004 99020 0.9948 exp 1 0
## 2 1e+05 0.06 1.062 1.062 94020 0.9993 exp 1 0
## 3 1e+05 0.11 1.119 1.118 89099 0.9928 exp 1 0
## 4 1e+05 0.16 1.191 1.196 83981 1.0024 exp 1 0
## 5 1e+05 0.21 1.255 1.253 79220 0.9905 exp 1 0
## 6 1e+05 0.26 1.351 1.353 73992 1.0008 exp 1 0
## label
## 1 exp(1,0)
## 2 exp(1,0)
## 3 exp(1,0)
## 4 exp(1,0)
## 5 exp(1,0)
## 6 exp(1,0)
Recommençons notre graphe
ggplot(df, aes(x = lambda, y = response, color = label)) + geom_line(aes(linetype = input_type)) +
geom_point() + geom_errorbar(width = 0.02, aes(x = lambda, y = response,
ymin = response - 2 * sd/sqrt(sys_vide), ymax = response + 2 * sd/sqrt(sys_vide))) +
geom_vline(xintercept = 1) + geom_hline(yintercept = 1) + theme_bw() + scale_color_brewer(palette = "Dark2") +
ylim(0, 6) + 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")
## Warning: Removed 12 rows containing missing values (geom_path).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
Bon, c'est bien plus propre maintenant. La loi déterministe est celle avec les temps de réponse les plus faibles et la loi gamma avec un paramètre de forme de 0.2 semble exploser bien avant que le taux d'arrivée atteigne 1. Après, difficile de dire si elle explose vraiment (i.e., si le système est vraiment instable et non définis) ou si le temps de réponse est juste très très grand.
La loi exponentielle est un peu entre ces deux extrêmes. Dans tous les cas, la variance du temps de service a un effet important sur le temps de réponse moyen. Essayons d'illustrer ceci en colorant les courbes en fonction de la variance du temps de service.
ggplot(df, aes(x = lambda, y = response, shape = label, colour = sd_S)) + geom_line() +
geom_point(size = 4) + scale_colour_continuous(low = "lightblue", high = "red") +
geom_errorbar(width = 0.03, aes(x = lambda, y = response, ymin = response -
2 * sd/sqrt(sys_vide), ymax = response + 2 * sd/sqrt(sys_vide))) + geom_vline(xintercept = 1) +
geom_hline(yintercept = 1) + theme_bw() + ylim(0, 6) + 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")
## Warning: Removed 12 rows containing missing values (geom_path).
## Warning: Removed 12 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
Parfait, on voit que la variance du temps de service a un impact fort sur le temps de réponse, et ce même avec des débits d'entrée faibles! D'autre part, certaines courbes sont quasiment superposées alors qu'elles correspondent à des lois différentes. En revanche elles correspondent à des distributions ayant à peu près la même variance… Il y a certainement un théorème qui dit que le temps de réponse d'une file FIFO M/G/1 ne dépend que du débit d'arrivée \( \lambda \), du débit de service \( \mu \) et de la variance du service. :)
D'autre part, quand on voit le temps que ça prend à calculer à cause des fortes variances, on peut se dire que ça vaudrait le coup d'appliquer les techniques de réduction de variance qu'on a vu précédemment! ;)