4.7 Análisis de tiempos de ocupación

En esta sección, nos concentramos en el análisis de los tiempos de ocupación de un estado determinado en un intervalo de tiempo finito \([0, T]\), es decir, la duración esperada de tiempo que el sistema pasa en ese estado. Como ocurre con las probabildiades de trasnsición presentamos un algortimo para obtener las cantidades de interés a partir de la matriz de tasas, y veremos como la simulación del proceso nos puede ayudar a dar respuesta a las mismas cuestiones.

Teorema 4.3 Si \(P(t)\) es la matriz de transiciones de probabilidad del proceso \(\{X(t), t \geq 0\}\) con espacio de estados \(S = \{1, 2,...,N\}\), se define la cantidad \(m_{ij}(T)\) como el tiempo de ocupación del estado \(j\) hasta el tiempo \(T\) partiendo del estado \(i\) como:

\[m_{ij}(T) = \int_0^T p_{ij}(t)dt, \quad 1 \leq i, j \leq N.\]

Cuando tenemos formas explicitas para cada uno de los elemntos de \(P(t)\) (este es el caso para la mayoría de los sistemas de colas de espera) este probelma se puede resolver teóricamente, pero en la mayoria de ocasiones es necesario un algoritmo de computación para aproximar estas cantidades. A continuación se presenta el algoritmo necesario, pero antes veamos la aproximación mediante series de la matriz \(M(T) = [m_{ij}(t)]\).

Teorema 4.4 Si \(Y\) es una variable aletoria Poisson de parámetro \(r*t\) entonces:

\[M(T) = \frac{1}{r} \sum_{k = 0}^{\infty} P(Y > k) \hat{P}^k, \quad T \geq 0.\]

Algoritmo para obtener \(M(T)\):

  1. Fijar \(R\), \(T\), y \(0 < \epsilon < 1\). Por defecto fijamos \(\epsilon = 0.00001\).
  2. Obtener r.
  3. Calcular \(\hat{P}\).
  4. Calcular \(A = \hat{P}\); \(k = 0\)
  5. Calcular \(yek = exp(-r*t)\), \(ygk = 1 - yek\), \(suma = ygk\)
  6. Calcular \(B = ygk * I\)
  7. Mientras que \(suma/r < T-\epsilon\), calcular:

\[k = k + 1; \quad yek = yek*(rT)/k; \quad ygk = ygk - yek\] \[B = B + ygk*A; \quad A = A\hat{P}; \quad suma = suma + ygk\] Al finalizar la matriz \(B/r\) es una aproximación de \(M(T)\) con error inferior a \(\epsilon\).

En realidad en nuestro algoritmo añadiremos un parámetro extra para indicar como se debe calcular el valor de \(r\), bien como el máximo los elementos por filas de \(R\) o como la suma de todos ellos.

tiempos.ocupacion<- function(Rmat, Ts, cal)
{
  # Algortimo de uniformización para obtener M(T)
  ################################################
  
  # Parámetros de la función
  # Rmat: matriz de tasas
  # ts: instante de tiempo
  # epsilon: error en la aproximación
  # cal: forms de obtener r con dos valores 1 = máximo, 2 = suma
  epsilon <- 1e-05
  # Paso 2. Calculo de r
  ris <- apply(Rmat, 1, sum)
  ifelse(cal == 1, rlimit <- max(ris), rlimit <- sum(Rmat))
  # Paso 3. Calculo de hat(P)
  hatP <- Rmat/rlimit
  diag(hatP) <- 1 - ris/rlimit
  # Paso 4. 
  k <- 0
  A <- hatP
  # Paso 5.
  yek <- exp(-1*rlimit*Ts)
  ygk <- 1 - yek
  suma <- ygk
  # Paso 6.
  B <- ygk*diag(nrow(Rmat))
  # Bucle simulación
  cota <- Ts- epsilon
  while(suma/rlimit < cota)
  {
    k <- k + 1
    yek <- yek*(rlimit*Ts)/k
    ygk <- ygk - yek
    B <- B + ygk*A
    A <- A%*%hatP
    suma <- suma + ygk
  }
  return(round(B/rlimit, 4))
}

Ejemplo 4.4 Retomamdo el sistema sobre el tiempo de vida de una máquina descrito en el ejemplo ??, supongamos que el tiempo esperado hasta que falla una máquina son 10 días, mientras que el tiempo esperado de reparación es de 1 día. Si la máquina funciona el primer día de enero ¿cuál es el tiempo total esperado de funcionamiento de la máquina al finalizar el mes de enero? Dado que el proceso sólo tiene espacio de estados \(S = \{0, 1\}\), la cantidad de interés es \(m_{11}(31)\). Con \(\lambda = 1\) y \(\mu = 0.1\), la matriz de tasas del proceso viene dada por:

nestados <- 2
lambda <- 1
mu <- 0.1
R <- matrix(nrow = nestados, ncol = nestados, data = 0)

R[1,2] <- lambda 
R[2,1] <- mu

Obtenemos ahora la matriz \(M(31)\) correspondiente a la cantidad de interés, teniendo en cuenta que partimos de que la máquina está en marcha (primera fila de \(M\)) al inicio:

tiempos.ocupacion(R, 31, 1)
##        [,1]    [,2]
## [1,] 3.6446 27.3554
## [2,] 2.7355 28.2645

Por tanto, el tiempo esperado de funcionamiento es de 28.26 días, mientras que el tiempo de reparación es de \(2.74\) días. Utilizando el simulador del proceso podemos ver que el resultado obtenido es similar.

# Réplicas del proceso
replicas <- 2500
envs <- lapply(1:replicas, function(i) {
    maquina <- sistema.1m(31, 1, 1/10)
})
# almacenamos análisis de llegadas del sistema
simarrivals <- as_tibble(get_mon_arrivals(envs))
# Almacenamos el estado final de la cola en el último instante del sistema
simarrivals %>%
  mutate(tOFF = end_time - start_time) %>%
  group_by(replication) %>%
  summarise(totalOFF = sum(tOFF)) %>%
  ungroup() %>%
  summarise(mOFF = mean(totalOFF), mON = 31 - mOFF)
## # A tibble: 1 × 2
##    mOFF   mON
##   <dbl> <dbl>
## 1  2.85  28.2

Para practicar este apartado puedes resolver los ejercicios B-5 a B-8 de la colección al final de la unidad.