1. Introducción

La teoría de colas estudia sistemas con llegadas aleatorias de usuarios (clientes, trabajos) a un servidor (o varios), donde los clientes pueden esperar si el servidor está ocupado. Notación estándar de Kendall: A/S/c (llegadas A, servicio S, c servidores), con A y S usualmente de la forma M (Poisson/Exponencial), D (determinístico), E_k (Erlang), G (general).

En esta nota abordamos: procesos nacimiento–muerte (BD), M/M/1, M/M/c, M/G/1 (fórmula P–K), Leyes de Little, PASTA, y métricas (\(L, L_q, W, W_q\)), más simulación en R.

2. Marco BD y estacionariedad

Sea \(\{N(t)\}_{t\ge 0}\) el número de clientes en sistema. En muchos modelos clásicos, \(N(t)\) es una CMTC de nacimiento–muerte con tasas: - Nacimiento (llegada): \(\lambda_n\) desde estado \(n\) a \(n+1\), - Muerte (servicio): \(\mu_n\) desde estado \(n\) a \(n-1\).

La matriz generadora \(Q\) cumple \(q_{n,n+1}=\lambda_n\), \(q_{n,n-1}=\mu_n\), \(q_{n,n}=-(\lambda_n+\mu_n)\).
La ecuación de balance detallado (suficiente para estacionariedad en BD irreducible) es \[ \pi_n \lambda_n = \pi_{n+1}\mu_{n+1}, \] de donde \[ \pi_n = \pi_0 \prod_{k=0}^{n-1}\frac{\lambda_k}{\mu_{k+1}},\quad \text{con } \sum_{n\ge 0}\pi_n=1. \] Cuando \(\sum_{n\ge 0}\prod_{k=0}^{n-1}\frac{\lambda_k}{\mu_{k+1}}<\infty\) existe distribución estacionaria.

3. Modelo M/M/1

Supuestos: Llegadas Poisson tasa \(\lambda>0\) (PASTA), tiempos de servicio Exp\((\mu)\) i.i.d., un servidor, disciplina FCFS. Entonces \(N(t)\) es BD con \[ \lambda_n=\lambda,\qquad \mu_n=\mu \ (n\ge 1),\quad \mu_0=0. \] Sea \(\rho=\lambda/\mu\) (utilización). Condición de estabilidad: \(\rho<1\).

3.1 Distribución estacionaria (demostración breve)

Balance detallado: \(\pi_n\lambda=\pi_{n+1}\mu\Rightarrow \pi_{n+1}=\rho\,\pi_n\).
Por inducción: \(\pi_n=\rho^n\pi_0\). Normalizando \(\sum_{n\ge 0}\pi_n=\pi_0\sum_{n\ge 0}\rho^n=1\) da \(\pi_0=1-\rho\), por tanto \[ \boxed{\ \pi_n=(1-\rho)\rho^n,\quad n\ge 0,\ \rho<1\ }. \]

3.2 Métricas clásicas

\[ L=\mathbb E[N]=\sum_{n\ge 0} n\pi_n=\frac{\rho}{1-\rho},\qquad L_q=\mathbb E[N_q]=\frac{\rho^2}{1-\rho}. \] Por Little (\(L=\lambda W\), \(L_q=\lambda W_q\)): \[ W=\frac{1}{\mu-\lambda},\qquad W_q=\frac{\lambda}{\mu(\mu-\lambda)}=\frac{\rho}{\mu-\lambda}. \]

Esbozo de prueba Little: En estado estable, en régimen de largo plazo, el número promedio en sistema \(L\) es flujo (tasa efectiva \(\lambda\)) por tiempo medio en sistema \(W\). Esto se formaliza con identidades de tiempo de permanencia (integrales de ocupación) y límites ergódicos.

4. Modelo M/M/c

Supuestos: \(c\) servidores en paralelo, llegadas Poisson \(\lambda\), servicios Exp\((\mu)\) i.i.d. Si hay \(n\ge c\) en el sistema, la tasa total de servicio es \(c\mu\); si \(n<c\), es \(n\mu\). Así, \[ \lambda_n=\lambda,\quad \mu_n=\min(n,c)\mu. \] Defina \(\rho=\frac{\lambda}{c\mu}\) (utilización por servidor). Estabilidad: \(\rho<1\).

La distribución estacionaria (vía BD) resulta: \[ \pi_n=\begin{cases} \displaystyle \pi_0\frac{(\lambda/\mu)^n}{n!}, & 0\le n<c,\\[6pt] \displaystyle \pi_0\frac{(\lambda/\mu)^n}{c!\,c^{\,n-c}}, & n\ge c, \end{cases} \] con \[ \pi_0=\left[\sum_{n=0}^{c-1}\frac{(\lambda/\mu)^n}{n!}+\frac{(\lambda/\mu)^c}{c!}\frac{1}{1-\rho}\right]^{-1}. \]

La prob. de espera (cola no vacía) es la Erlang C: \[ P_{\text{wait}}=\Pr\{N\ge c\}= \frac{\frac{(\lambda/\mu)^c}{c!\,}\frac{1}{1-\rho}}{\sum_{n=0}^{c-1}\frac{(\lambda/\mu)^n}{n!}+\frac{(\lambda/\mu)^c}{c!\,}\frac{1}{1-\rho}}. \] Luego \[ L_q=\frac{P_{\text{wait}}\ \rho}{1-\rho},\qquad W_q=\frac{L_q}{\lambda},\qquad W=W_q+\frac{1}{\mu},\qquad L=\lambda W. \]

5. Modelo M/G/1 y fórmula de Pollaczek–Khinchine (P–K)

Ahora un servidor, llegadas Poisson \(\lambda\), tiempos de servicio generales \(S\) con media \(\mathbb E[S]=1/\mu\) y segunda momento finito \(\mathbb E[S^2]\). Defina \(\rho=\lambda \mathbb E[S]\) (estabilidad: \(\rho<1\)).

La fórmula P–K (para FCFS) da: \[ \boxed{\ L_q=\frac{\lambda^2 \,\mathbb E[S^2]}{2(1-\rho)}\ ,\quad W_q=\frac{L_q}{\lambda}=\frac{\lambda\,\mathbb E[S^2]}{2(1-\rho)},\quad W=W_q+\mathbb E[S],\ L=\lambda W\ .} \]

Esbozo de demostración: Usar el árbol de residuos (residual life) y el método de tiempos de ocupación: al llegar un cliente, por PASTA, ve el tiempo residual de servicio del cliente en curso con esperanza \(\mathbb E[S^2]/(2\mathbb E[S])\), más el trabajo delante (suma de servicios de los que esperan); linealizando esperanzas y usando estabilidad se obtiene P–K.

6. PASTA y Little (demostraciones breves)

PASTA (Poisson Arrivals See Time Averages): si las llegadas son Poisson y el proceso subyacente es independiente de la futura llegada, entonces lo que ven los arribos coincide con promedios temporales estacionarios (prueba vía intensidades y compensadores; los instantes de llegada muestrean el proceso sin sesgo temporal).
Little: para cualquier sistema estable, sin supuestos de distribución, \[ \boxed{L=\lambda W,\qquad L_q=\lambda W_q.} \] Prueba mediante identidad de ocupación: el área bajo \(N(t)\) en un horizonte \([0,T]\) equivale a la suma de los tiempos en sistema de todos los atendidos, y al tomar límites \(T\to\infty\) se obtienen identidades de tasa × tiempo medio.

7. Simulación en R: M/M/1 y M/M/c

7.1 Simulación M/M/1 (Gillespie/next-event)

set.seed(123)

simulate_mm1 <- function(lambda, mu, T_end, x0 = 0){
  t <- 0; n <- x0
  times <- c(0); levels <- c(n)
  while(t < T_end){
    rate_arr <- lambda
    rate_dep <- ifelse(n > 0, mu, 0)
    rate_tot <- rate_arr + rate_dep
    if(rate_tot == 0) break
    dt <- rexp(1, rate_tot)
    t_next <- t + dt
    if(t_next > T_end){ t <- T_end; times <- c(times, t); levels <- c(levels, n); break }
    # evento
    if(runif(1) < rate_arr/rate_tot){ n <- n + 1 } else { n <- n - 1 }
    t <- t_next
    times <- c(times, t); levels <- c(levels, n)
  }
  data.frame(time = times, N = levels)
}

# ejemplo
lambda <- 0.7; mu <- 1
rho <- lambda/mu
path <- simulate_mm1(lambda, mu, T_end = 200, x0 = 0)

# Estimación empírica de L = promedio temporal de N(t)
dt <- diff(path$time); midN <- path$N[-1]
L_emp <- sum(dt * midN) / sum(dt)
L_theo <- rho/(1-rho)
c(L_emp = round(L_emp,3), L_theo = round(L_theo,3))

##  L_emp L_theo 
##  3.029  2.333

simulate_mm1_fcfs_waits <- function(lambda, mu, n_customers = 1e4, warmup = 1000){
  arr <- cumsum(rexp(n_customers, lambda))    # tiempos de llegada
  serv <- rexp(n_customers, mu)               # tiempos de servicio
  start <- numeric(n_customers)               # inicio de servicio
  finish <- numeric(n_customers)              # fin de servicio
  
  start[1] <- arr[1]
  finish[1] <- start[1] + serv[1]
  for(i in 2:n_customers){
    start[i] <- max(arr[i], finish[i-1])
    finish[i] <- start[i] + serv[i]
  }
  wait_q <- pmax(0, start - arr)              # Wq
  wait   <- finish - arr                      # W
  idx <- (warmup+1):n_customers               # descartar calentamiento
  c(Wq = mean(wait_q[idx]), W = mean(wait[idx]))
}

set.seed(321)
res <- simulate_mm1_fcfs_waits(lambda, mu)
W_theo  <- 1/(mu - lambda)
Wq_theo <- lambda/(mu*(mu - lambda))
c(sim_Wq = round(res["Wq"],3), theo_Wq = round(Wq_theo,3),
  sim_W  = round(res["W"],3),  theo_W  = round(W_theo,3))

## sim_Wq.Wq   theo_Wq   sim_W.W    theo_W 
##     2.391     2.333     3.385     3.333

simulate_mm1_fcfs_waits <- function(lambda, mu, n_customers = 1e4, warmup = 1000){
  arr <- cumsum(rexp(n_customers, lambda))    # tiempos de llegada
  serv <- rexp(n_customers, mu)               # tiempos de servicio
  start <- numeric(n_customers)               # inicio de servicio
  finish <- numeric(n_customers)              # fin de servicio
  
  start[1] <- arr[1]
  finish[1] <- start[1] + serv[1]
  for(i in 2:n_customers){
    start[i] <- max(arr[i], finish[i-1])
    finish[i] <- start[i] + serv[i]
  }
  wait_q <- pmax(0, start - arr)              # Wq
  wait   <- finish - arr                      # W
  idx <- (warmup+1):n_customers               # descartar calentamiento
  c(Wq = mean(wait_q[idx]), W = mean(wait[idx]))
}

set.seed(321)
res <- simulate_mm1_fcfs_waits(lambda, mu)
W_theo  <- 1/(mu - lambda)
Wq_theo <- lambda/(mu*(mu - lambda))
c(sim_Wq = round(res["Wq"],3), theo_Wq = round(Wq_theo,3),
  sim_W  = round(res["W"],3),  theo_W  = round(W_theo,3))

## sim_Wq.Wq   theo_Wq   sim_W.W    theo_W 
##     2.391     2.333     3.385     3.333

erlangC <- function(lambda, mu, c){
  a <- lambda/mu
  rho <- lambda/(c*mu)
  s1 <- sum(sapply(0:(c-1), function(n) a^n/factorial(n)))
  s2 <- a^c/factorial(c) * 1/(1 - rho)
  P0 <- 1/(s1 + s2)
  Pw <- s2 * P0
  Lq <- Pw * rho/(1 - rho)
  Wq <- Lq/lambda
  W  <- Wq + 1/mu
  L  <- lambda * W
  list(P0=P0, Pw=Pw, Lq=Lq, Wq=Wq, W=W, L=L, rho=rho)
}

lambda <- 4; mu <- 3; c <- 2
out <- erlangC(lambda, mu, c)
lapply(out, function(x) round(x,4))

## $P0
## [1] 0.2
## 
## $Pw
## [1] 0.5333
## 
## $Lq
## [1] 1.0667
## 
## $Wq
## [1] 0.2667
## 
## $W
## [1] 0.6
## 
## $L
## [1] 2.4
## 
## $rho
## [1] 0.6667

set.seed(999)

simulate_mg1_fcfs <- function(lambda, r_shape, r_rate, n_customers = 3e4, warmup = 2e3){
  arr <- cumsum(rexp(n_customers, lambda))
  serv <- rgamma(n_customers, shape = r_shape, rate = r_rate) # E[S]=shape/rate
  start <- numeric(n_customers); finish <- numeric(n_customers)
  start[1] <- arr[1]; finish[1] <- start[1] + serv[1]
  for(i in 2:n_customers){
    start[i] <- max(arr[i], finish[i-1])
    finish[i] <- start[i] + serv[i]
  }
  Wq <- pmax(0, start - arr)
  W  <- finish - arr
  idx <- (warmup+1):n_customers
  c(Wq = mean(Wq[idx]), W = mean(W[idx]), ES = mean(serv), ES2 = mean(serv^2))
}

lambda <- 0.8
shape <- 3; rate <- 3   # E[S]=shape/rate=1,  Var=shape/rate^2=1/3
res <- simulate_mg1_fcfs(lambda, shape, rate)
rho <- lambda * res["ES"]
Wq_PK <- (lambda * res["ES2"]) / (2 * (1 - rho))
c(sim_Wq = round(res["Wq"],3), PK_Wq = round(Wq_PK,3),
  ES = round(res["ES"],3), rho = round(rho,3))

## sim_Wq.Wq PK_Wq.ES2     ES.ES    rho.ES 
##     2.584     2.699     1.003     0.802

Teoría de Colas

Jhonier Rangel

2025-09-17