Licenciado bajo Creative Commons BY-NC-SA 4.0.

1 Sistemas de colas

En este tutorial, continuaremos utilizando técnicas DES y el paquete simmer. Consulte el tutorial previo si necesita recordar la Introducción a DES con R.

library(simmer)

## Warning: package 'simmer' was built under R version 3.6.3

library(simmer.plot)

## Warning: package 'simmer.plot' was built under R version 3.6.3

1.1 Sistema M/M/1

Un sistema M/M/1 es aquel con llegadas exponenciales (M/M/1), un servidor (M/M/1) con tiempo de servicio exponencial (M/M/1) y cola infinita (M/M/1/\(\infty\) implícito).

Ejemplo 1: Un cajero al que llegan personas según un proceso de Poisson de tasa \(\lambda\), hacen cola en la calle y, cuando llega su turno, sacan dinero según otro proceso de Poisson de tasa \(\mu\).

Recordemos los parámetros básicos de este sistema:

\[\begin{aligned} \rho &= \frac{\lambda}{\mu} &&\equiv \mbox{Ocupación media del servidor} \\ N &= \frac{\rho}{1-\rho} &&\equiv \mbox{Número medio de usuarios en el sistema (cola + servidor)} \\ T &= \frac{N}{\lambda} &&\equiv \mbox{Tiempo medio de estancia en el sistema (cola + servidor) [teorema de Little]} \\ \end{aligned}\]

Siempre que \(\rho < 1\). Si esto no se cumple, quiere decir que en promedio hay más llegadas de las que el servidor puede absorber, por lo que la cola crecerá indefinidamente.

La simulación de un sistema M/M/1 es inmediata partiendo del Ejemplo 1:

# Tomamos estos valores para lambda y mu
lambda <- 2
mu <- 4
rho <- lambda/mu # = 2/4

# Trayectoria que siguen las personas
m.queue <- trajectory() %>%
  seize("server", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("server", amount=1)

# Simulador en el que incluimos un cajero y una cola infinita, junto con un generador usando
## la trayectoria previamente definida
mm1.env <- simmer() %>%
  add_resource("server", capacity=1, queue_size=Inf) %>%
  add_generator("arrival", m.queue, function() rexp(100, lambda)) %>%
  run(4000/lambda)

# Simulamos y almacenamos los resultados
mm1.df.res <- get_mon_resources(mm1.env)
mm1.df.arr <- get_mon_arrivals(mm1.env)

# Valor teórico del número medio de usuarios en el sistema
mm1.N <- rho/(1-rho)

# Evolución del número medio de usuarios en el sistema + Valor teórico
plot(mm1.df.res, metric="usage", "server", items="system") + 
  geom_hline(yintercept = mm1.N)

De forma experimental, obtenemos el tiempo de estancia de cada usuario. A partir de estos valores, podemos obtener el tiempo medio, que coincide (aprox.) con el valor teórico:

mm1.t_system <- mm1.df.arr$end_time - mm1.df.arr$start_time
mean(mm1.t_system) ; mm1.N/lambda

## [1] 0.4830735

## [1] 0.5

La inversa de la diferencia media entre llegadas es la tasa efectiva, que en este sistema coincide (aprox.) con la lambda real porque no hay rechazo:

mm1.df.arr.finished <- subset(mm1.df.arr, finished == TRUE)

# Tasa efectiva: 1 / media de la diferencia de tiempo entre llegadas
1/mean(diff(mm1.df.arr.finished$start_time)) ; lambda

## [1] 1.952691

## [1] 2

# Tasa de rechazo: 1 - numero de usuarios aceptados por el sistema / total de usuarios
1 - nrow(mm1.df.arr.finished) / nrow(mm1.df.arr)

## [1] 0

Además, en un M/M/1 se cumple que la distribución de tiempos de estancia es a su vez una variable exponencial de media \(T\):

qqplot(mm1.t_system, rexp(1000, lambda/mm1.N))
abline(0, 1, lty=2, col="red")

1.2 Sistemas M/M/1/k

El sistema M/M/1/K es muy similar al M/M/1, pero el número máximo de usuarios en el sistema está limitado por K. Por lotanto, la longitud máxima de la cola esK−1, dado que el recurso puede albergar a un usuario. Este sistema resulta apropiado paramodelar situaciones de capacidad finita (una línea de transmission con un buffer de tamaño limitado, una sala de espera con aforo limiatdo, …). Partiendo del Ejemplo 1, definimos un sistema con 1 recurso y un máximo de 2 usuarios en el sistema:

lambda <- 2
mu <- 4

# total de usuarios en el sistema = usuarios en cola + usuarios en los recursos
k <- 2

# Número de recursos
m12_capacity = 1

# Capacidad de la cola
m12_queue_size = 1

# Usamos la trayectoria definida anteriormente 
m12.queue <- trajectory() %>%
  seize("server", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("server", amount=1)

# Simulador
mm12.env <- simmer() %>%
  add_resource("server", capacity=m12_capacity, queue_size=m12_queue_size) %>%
  add_generator("arrival", m12.queue, function() rexp(200, lambda)) %>%
  run(4000/lambda)

# Lanzamos las simulaciones
mm12.df.res <- get_mon_resources(mm12.env)
mm12.df.arr <- get_mon_arrivals(mm12.env)

1.2.1 Probabilidades

# Sistema vacío
## Teórica
I <- lambda/mu
n <- 0
pn <- ((1-I)*(I^n))/(1-I^(k+1))
pn

## [1] 0.5714286

## Simmer
deltas <- diff(mm12.df.res$time)
busy <- which(mm12.df.res$system == 0)
t_busy <- deltas[busy]
sum(t_busy, na.rm=TRUE) / max(mm12.df.res$time)

## [1] 0.5815656

# Sistema lleno
## Teórica
I <- lambda/mu
n <- k
pk <- ((1-I)*(I^n))/(1-I^(k+1))
pk

## [1] 0.1428571

## Simmer
deltas <- diff(mm12.df.res$time)
busy <- which(mm12.df.res$system == 2)
t_busy <- deltas[busy]
pk_simmer <- sum(t_busy, na.rm=TRUE) / max(mm12.df.res$time)
pk_simmer

## [1] 0.1419972

## or
pk_simmer_2 <- sum(!mm12.df.arr$finished) / nrow(mm12.df.arr)
pk_simmer_2

## [1] 0.1425652

1.2.2 Tasa efectiva

# Teórica
lambda_effect <- lambda*(1-pk)
lambda_effect

## [1] 1.714286

# simmer
lambda_effect_simmer <- lambda * (1 - pk_simmer)
lambda_effect_simmer

## [1] 1.716006

## or
mm12.df.arr.finished <- subset(mm12.df.arr, finished == TRUE)
1/mean(diff(mm12.df.arr.finished$start_time))

## [1] 1.678551

1.2.3 Número medio de usuarios

# Teórico
N <- (I/(1-I)) - (k+1)*(I^(k+1))/(1-(I^(k+1)))
N

## [1] 0.5714286

# simmer
deltas_2 <- diff(mm12.df.res$time)
deltas_2 <- c(deltas_2, 0)
## Ponderamos el numero de usuarios en cada instante por el tiempo que hemos tenido ese número de usuarios
res <- mm12.df.res$system * deltas_2
## Hacemos la media del numero de usuarios ponderado
N_simmer <- sum(res, na.rm=TRUE) / max(mm12.df.res$time)
N_simmer

## [1] 0.5601368

1.2.4 Tiempo medio de estancia

# Teórico
T <- N/lambda_effect
T

## [1] 0.3333333

#simmer
mm12.t_system <- mm12.df.arr$end_time - mm12.df.arr$start_time
mean(mm12.t_system)

## [1] 0.2861962

1.3 Sistemas M/M/m/k

Un sistema M/M/c/k mantiene las llegadas y tiempos de servicio exponenciales, pero consta de más de un servidor en general y una cola de tamaño limitado, lo que nos permite modelar problemas más realistas. Por ejemplo, un router podrá disponer de varios procesadores para reenviar paquetes y las colas que almacenan paquetes serán necesariamente de tamaño finito.

Simulamos a continuación un sistema M/M/2/4 (2 servidores, 2 posiciones en cola) siguiendo con el Ejemplo 1. Obsérvese que la trayectoria no varía con respecto al caso M/M/1 y M/M/1/k.

lambda <- 2
mu <- 4
m <- 2
capacity <- 2
queue_size <- 2
k <- 4

# Usamos la trayectoria definida anteriormente 
m13.queue <- trajectory() %>%
  seize("server", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("server", amount=1)

mm13.env <- simmer() %>%
  add_resource("server", capacity=capacity, queue_size=queue_size) %>%
  add_generator("arrival", m13.queue, function() rexp(100, lambda)) %>%
  run(4000/lambda)

mm13.df.res <- get_mon_resources(mm13.env)
mm13.df.arr <- get_mon_arrivals(mm13.env)

De manera similar al apartado anterior podremos obtener probabilidades, tasa efectiva, número medio de usuarios en el sistema y tiempo medio de estancia en el sistema.

2 Ejemplo: Diseño de redes de comunicaciones

2.1 Planteamiento

Considere un operador de telefonía móvil que desea dar cobertura a un área de 10 km\(^2\). Para ello, deberá dividir el área en \(N\) parcelas, en cada una de las cuales colocará una estación base. Se dispone de dos tipos de estación base en función del número de canales (llamadas que puede cursar simultáneamente) y el coste:

Estación base	A	B
Número de canales	2	3
Coste por unidad (€)	10.000	15.000

cells

Se sabe que los usuarios están distribuidos uniformemente con una densidad de 10 usuarios por km\(^2\). Cada uno de ellos genera llamadas según un proceso de Poisson de tasa 2 llamadas por hora en la hora pico. Por tanto,

\[\lambda = 10 \cdot 10 \cdot 2 = 200\]

Las llamadas tienen una duración media de 3 minutos:

\[\mu = \frac{60}{3} = 20\]

Consideramos que, cuando una llamada no dispone de recursos para ser atendida, se rechaza, y que el usuario no reintenta llamar. En estas condiciones, una estación base se puede modelar mediante un sistema M/M/m/m, donde \(m\) es el número de canales disponibles y no hay cola de espera.

Se desea decidir mediante simulación qué tipo de estación base se va a desplegar (a) o (b). Si se impone una tasa de rechazo (llamadas que no llegan a ser atendidas) inferior al 5 %, se trata de hallar el número de estaciones base necesarias al mínimo coste.

2.2 Resolucion

En primer lugar vamos a definir una función que, dado un número de canales, nos diga el número de estaciones que necesitaríamos para cumplir la tasa de rechazo requerida. Para ello, simulamos una celda a la que llega una intensidad de tráfico igual a \(\lambda\) dividido por el número de estaciones desplegadas. A continuación, realizamos 10 búsquedas para cada tipo de estación considerando unas 1000 llamadas en cada simulación:

station_a_channels <- 2
station_a_cost <- 10000

station_b_channels <- 3
station_b_cost <- 15000

lambda <- 200
mu <- 20
# Consideramos 1000 llamadas en cada simulacion
calls <- 1000

m14.queue <- trajectory() %>%
  seize("server", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("server", amount=1)

find_n_stations <- function(channels, start=1) {
  
  stations <- start
  
  while (TRUE) {
    cell.arr <- simmer() %>%
      add_resource("server", capacity=channels, queue_size=0) %>%
      add_generator("arrival", m14.queue, function() rexp(100, lambda/stations)) %>%
      run(calls*stations/lambda) %>%
      get_mon_arrivals()
    
    # Tasa de rechazo
    rej_rate <- sum(!cell.arr$finished) / nrow(cell.arr)
    
    if (rej_rate < 0.05){
      break;
    } 
    
    stations <- stations + 1
  }
  return(stations)
}

c_a <- c()
c_b <- c()
for( i in c(1:10)){
  c_a <- c(c_a, find_n_stations(station_a_channels))
  c_b <- c(c_b, find_n_stations(station_b_channels))
}
avg_a_stations <- mean(c_a)
avg_b_stations <- mean(c_b)

avg_a_stations

## [1] 25.6

avg_b_stations

## [1] 11.1

Podemos comprobar que el número de estaciones medio es superior a 20 para el caso (a) y superior a 10 para el caso (b). Por tanto, podemos utilizar estos valores en el parámetro start para optimizar la búsqueda.

c_a <- c()
c_b <- c()
for( i in c(1:10)){
  c_a <- c(c_a, find_n_stations(station_a_channels, start = 20))
  c_b <- c(c_b, find_n_stations(station_b_channels, start = 10))
}
avg_a_accurate <- mean(c_a)
avg_b_accurate <- mean(c_b)
avg_a_accurate; avg_b_accurate

## [1] 25.5

## [1] 11.4

cost_a <- avg_a_accurate * 10000
cost_b <- avg_b_accurate * 15000

cost_a; cost_b

## [1] 255000

## [1] 171000

El resultado, de media, es que se requieren más de 25 estaciones en el caso a) y más de 11 estaciones en el caso b), que además resulta el más barato.

3 Resolución de problemas con R

3.1 Problema 7.7

Suponga una barra de comida rápida en la que los clientes se sientan en uno de los dos taburetes disponibles a comer. Los clientes llegan según un proceso de Poisson de tasa 10 clientes/hora, si no hay taburete disponible se marchan, y el tiempo que pasan en la barra se puede modelar con una variable aleatoria exponencial de media 12 minutos. En estas condiciones, cada cliente paga en media 30€. El dueño del local se plantea ampliar la barra con otro taburete, pero cree que los clientes pagarían 5€ menos cada uno. Calcule:

Probabilidad de que la barra esté vacía en las dos situaciones

¿Ganaría más dinero ampliando la barra?

#M/M/Resources/capacity

lambda <- 10
mu <- 5
I <- lambda/mu

m77.queue <- trajectory() %>%
  seize("taburete", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("taburete", amount=1)

############################
############ a) ############
############################

## 2 taburetes

mm77_2.df.res <- lapply(1:50, function(i) {
  simmer() %>%
    add_resource("taburete", capacity=2, queue_size=0) %>%
    add_generator("arrival", m77.queue, function() rexp(100, lambda)) %>%
    run(2000/lambda)
}) %>% get_mon_resources()

## Probabilidad de que el sistema esté vacío (teorico = 0.2)
res <- c()
for(i in c(1:50)){
  temp <- mm77_2.df.res[which(mm77_2.df.res$replication == i), ]
  
  deltas <- diff(temp$time)
  busy <- which(temp$system == 0)
  t_busy <- deltas[busy]
  res <- c(res, sum(t_busy, na.rm=TRUE) / max(temp$time))
}
mean(res)


## 3 taburetes
mm77_3.df.res <- lapply(1:50, function(i) {
  simmer() %>%
    add_resource("taburete", capacity=3, queue_size=0) %>%
    add_generator("arrival", m77.queue, function() rexp(100, lambda)) %>%
    run(2000/lambda)
}) %>% get_mon_resources()

## Probabilidad de que el sistema esté vacío (teorico = 0.157)
res <- c()
for(i in c(1:50)){
  temp <- mm77_3.df.res[which(mm77_3.df.res$replication == i), ]
  deltas <- diff(temp$time)
  busy <- which(temp$system == 0)
  t_busy <- deltas[busy]
  res <- c(res, sum(t_busy, na.rm=TRUE) / max(temp$time))
}
mean(res)

############################
############ b) ############
############################

## 2 taburetes

mm2.df.arr <- lapply(1:50, function(i) {
  simmer() %>%
    add_resource("taburete", capacity=2, queue_size=0) %>%
    add_generator("arrival", m77.queue, function() rexp(100, lambda)) %>%
    run(2000/lambda)
}) %>% get_mon_arrivals()

## tasa de rechazo es el número de clientes que no ha podido entrar entre el número de clientes que han llegado
rej_rate_2 <- sum(!mm2.df.arr$finished) / nrow(mm2.df.arr)

### Earnings =  tasa efectiva de clientes * precio
lambda * (1 - rej_rate_2) * 30

## 3 taburetes
mm3.df.arr <- lapply(1:50, function(i) {
  simmer() %>%
    add_resource("taburete", capacity=3, queue_size=0) %>%
    add_generator("arrival", m77.queue, function() rexp(100, lambda)) %>%
    run(2000/lambda)
}) %>% get_mon_arrivals()

## tasa de rechazo es el número de clientes que no ha podido entrar entre el número de clientes que han llegado
rej_rate_3 <- sum(!mm3.df.arr$finished) / nrow(mm3.df.arr)

### Earnings =  tasa efectiva de clientes * precio
lambda * (1 - rej_rate_3) * 25

3.2 Problema 7.8

Una pequeña central de comunicaciones da servicio a un municipio de 30.000 habitantes. Se estima que cada habitante realiza, en media, una llamada cada 30 días, y que la duración de estas llamadas se puede modelar con una v.a. exponencial de media 4 minutos y 30 segundos. Si la central dispone de 6 circuitos, por lo que puede cursar como máximo 6 llamadas a la vez, calcule:

Probabilidad de que una llamada sea rechazada.

Número medio de circuitos ocupados.

#M/M/Resources/capacity

lambda <- 125/3
mu <- 40/3
rho <- lambda/mu

m78.queue <- trajectory() %>%
  seize("circuito", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("circuito", amount=1)

############################
############ a) ############
############################

mm78.df.res <- lapply(1:50, function(i) {
  simmer() %>%
    add_resource("circuito", capacity=6, queue_size=0) %>%
    add_generator("arrival", m78.queue, function() rexp(200, lambda)) %>%
    run(2000/lambda)
}) %>% get_mon_resources()

### Probabilidad de que el sistema esté lleno (teórico = 0.059)
res <- c()
for(i in c(1:50)){
  temp <- mm78.df.res[which(mm78.df.res$replication == i), ]
  deltas <- diff(temp$time)
  busy <- which(temp$system == 6)
  t_busy <- deltas[busy]
  res <- c(res, sum(t_busy, na.rm=TRUE) / max(temp$time))
}
mean(res)

############################
############ b) ############
############################

## Por cada iteración, calculamos el número medio de circuitos ocupados (Teórico 2.94)
results <- c()
for(i in c(1:50)){
  temp <- mm78.df.res[which(mm78.df.res$replication == i), ]
  
  deltas_78 <- diff(temp$time)
  deltas_78 <- c(deltas_78, 0)
  ## Ponderamos el numero de usuarios en cada instante por el tiempo que hemos tenido ese número de usuarios
  res <- temp$system * deltas_78
  ## Hacemos la media del numero de usuarios ponderado
  results <- c(results, (sum(res, na.rm=TRUE) / max(temp$time)))
}
mean(results)

Teoría de Redes - Práctica 4

Iñaki Úcar, Pablo Serrano, José Alberto Hernández, Ginés García Avilés

22 de abril de 2020