Descripción

El ejemplo consiste en un sistema de colas con un único punto de servicio al que llegan clientes con tiempos entre llegadas dados por una distribución exponencial con media de 3 por hora. Los tiempos de servicio siguen una distribución exponencial con una duración media de 15 minutos.

Las medidas de rendimiento serían:

Una buena práctica de programación es importar todos los módulos y funciones necesarios al principio del programa.

library(knitr)
library(lubridate)

Inicialización

Las variables de estado son:

set.seed(12345678, "Mersenne-Twister")
lambda <- 3
mu <- 4
t <- 0
ocupado <- FALSE
cola <- 0
llegada <- rexp(1, lambda / 60)
prox_llegada <- t + llegada
prox_salida <- 9999
tsim <- 480

Los vectores tm, oc, y cl, almacenan el estado del sistema cada vez que ocurre una llegada o una salida, es decir, un evento.

Son respectivamente: - tm, el tiempo del evento en minutos. - oc, si el sistema está ocupado o no, True or False. - cl, es el número de clientes en la cola.

tm <- numeric()
oc <- logical()
cl <- integer()

clientes es una data frame donde se almacena la trayectoria de cada cliente en cuatro columnas

clientes <- data.frame("t_ll"=numeric(),
                       "t_sa"=numeric(),
                       "t_se"=numeric(),
                       "d_se"=numeric())

Bucle de la simulación

while (t < tsim){
  tm <- c(tm, t)
  oc <- c(oc, ocupado)
  cl <- c(cl, cola)
  if (prox_llegada < prox_salida){
    t <- prox_llegada
    llegada <- rexp(1, lambda / 60)
    prox_llegada <- t + llegada
    cliente <- data.frame("t_ll"=t, "t_sa"=NA, "t_se"=NA, "d_se"=NA)
    clientes <- rbind(clientes, cliente)
    if (ocupado == TRUE){
      cola <- cola + 1
    } else {
      ocupado <- TRUE
      servicio <- rexp(1, mu / 60)
      prox_salida <- t + servicio
      clientes[nrow(clientes), "t_sa"] = prox_salida
      clientes[nrow(clientes), "t_se"] = t
      clientes[nrow(clientes), "d_se"] = servicio
    }

  } else {
    t <- prox_salida
    if (cola > 0){
      cola <- cola - 1
      servicio <- rexp(1, mu / 60)
      prox_salida <- t + servicio
      clientes[nrow(clientes) - cola, "t_sa"] = prox_salida
      clientes[nrow(clientes) - cola, "t_se"] = t # la hora que lo empiezan a atender
      clientes[nrow(clientes) - cola, "d_se"] = servicio
    } else {
      ocupado <- FALSE
      prox_salida <- 9999
    }
  }
}

Reportes

Generar la tabla de clientes

cl_table <- data.frame("n_cli"=integer(),
                       "llega"=character(),
                       "sale"=character(),
                       "d_ocio"=numeric(),
                       "d_cola"=numeric(),
                       "d_serv"=numeric(),
                       "d_sis"=numeric())

t_ini = strptime("09:00 am", format="%I:%M %p")

time_str <- function(t_ini, t_ev){
  hr <- hours(t_ev %/% 60)
  min <- minutes(round(t_ev %% 60))
  sec <- seconds(round((t_ev * 60) %% 60))
  as.character(t_ini + hr + min + sec, format = "%H:%M:%S")
}

anterior <- 0
for (i in 1:nrow(clientes)){
  cliente <- clientes[i,]
  if (all(!is.na(cliente))){
    t_ll <- time_str(t_ini, cliente$t_ll)
    t_sa <- time_str(t_ini, cliente$t_sa)
    d_ocio <- max(0, cliente$t_ll - anterior)
    d_cola <- cliente$t_se - cliente$t_ll
    d_serv <- cliente$d_se
    d_sis <- cliente$t_sa - cliente$t_ll 
  }
  anterior <- cliente$t_sa
  row <- data.frame("n_cli"=i, "llega"=t_ll, "sale"=t_sa,
                    "d_ocio"=d_ocio, "d_cola"=d_cola, "d_serv"=d_serv, "d_sis"=d_sis)
  cl_table <- rbind(cl_table, row)
}
kable(cl_table, digits = 2, caption = "Traza de los clientes")
Traza de los clientes
n_cli llega sale d_ocio d_cola d_serv d_sis
1 09:04:00 09:23:19 4.00 0.00 19.32 19.32
2 09:04:07 09:25:43 0.00 19.20 1.40 20.60
3 09:05:24 09:35:50 0.00 19.32 10.12 29.44
4 09:36:15 09:44:32 1.41 0.00 7.29 7.29
5 09:49:13 10:03:37 5.68 0.00 13.39 13.39
6 09:52:43 10:23:45 0.00 10.90 20.14 31.04
7 10:08:27 10:42:32 0.00 14.31 18.78 33.09
8 10:35:42 10:55:09 0.00 6.83 13.62 20.45
9 11:00:28 11:01:08 5.32 0.00 0.66 0.66
10 11:08:15 11:11:46 7.12 0.00 2.52 2.52
11 11:13:31 11:33:16 1.74 0.00 20.76 20.76
12 11:16:19 11:58:14 0.00 16.95 24.97 41.92
13 11:20:14 12:19:11 0.00 38.01 20.95 58.96
14 12:36:03 12:41:35 16.87 0.00 4.53 4.53
15 12:47:36 12:52:40 6.02 0.00 5.07 5.07
16 12:57:43 13:00:47 5.05 0.00 3.06 3.06
17 13:08:30 13:26:09 8.72 0.00 17.65 17.65
18 13:21:15 13:39:04 0.00 4.90 12.92 17.81
19 15:35:20 16:14:21 116.28 0.00 39.01 39.01
20 16:03:14 16:15:17 0.00 11.11 0.92 12.04
21 16:04:24 16:20:49 0.00 10.88 4.54 15.42
22 16:04:28 16:21:54 0.00 15.35 1.08 16.44
23 16:21:32 16:53:17 0.00 0.36 32.39 32.74
24 16:36:05 17:03:52 0.00 17.19 9.59 26.78
25 16:36:05 17:03:52 0.00 17.19 9.59 26.78
26 16:36:05 17:03:52 0.00 17.19 9.59 26.78

Gráficos

Número de clientes en el sistema.

plot(tm, cl + oc, col="green", type="s", main="Número de Clientes en el sistema",
     xlab = "tiempo (min)", ylab = "# Clientes")

Número de clientes en la cola.

plot(tm, cl, col="red", type="s", main="Número de Clientes en la Cola",
     xlab = "tiempo (min)", ylab = "# Clientes")

Medidas de Rendimiento

Número promedio de clientes en el sistema

prod <- numeric(length(tm) - 1)
for (i in 1:length(tm)-1){
  prod[i] <- (tm[i+1] - tm[i]) * (cl[i] + oc[i])
}
L <- sum(prod) / tm[length(tm)]

\(L =\) 1.0372773.

Número promedio de clientes en la cola

prod <- numeric(length(tm) - 1)
for (i in 1:length(tm)-1){
  prod[i] <- (tm[i+1] - tm[i]) * (cl[i])
}
Lq <- sum(prod) / tm[length(tm)]

\(L_q =\) 0.4138171.

Fracción estimada de tiempo ocioso

prod <- numeric(length(tm) - 1)
for (i in 1:length(tm)-1){
  prod[i] <- (tm[i+1] - tm[i]) * (!oc[i])
}
P0 <- sum(prod) / tm[length(tm)]

\(P_0 =\) 0.3765399.

Ocupación promedio estimada

prod <- numeric(length(tm) - 1)
for (i in 1:length(tm)-1){
  prod[i] <- (tm[i+1] - tm[i]) * (oc[i])
}
Rho <- sum(prod) / tm[length(tm)]

\(\rho =\) 0.6234601.

Tiempo de espera promedio en el sistema

W <- sum(cl_table$d_sis) / nrow(cl_table)

\(W =\) 20.9047185

Tiempo de espera promedio en la cola

Wq <- sum(cl_table$d_cola) / nrow(cl_table)

\(W_q =\) 8.4493236