En este tutorial se abordará el tema de análisis de tiempos a través de un ejemplo. Para esto, utilizaremos la librería markovchain.
Vamos a resolver el problema del archivo “Complementaria 9 (Q).pdf” que se encuentra en SICUA+.
En primer lugar, se modelará la situación como una cadena de Markov de tiempo continuo. Se definen tres variables de estado:
\[\begin{eqnarray*} &X(t) = \text{ el número de unidades en la zona de prensas de presión en el tiempo }t \\ &Y(t) = \text{ el número de unidades en la zona de recubrimiento en el tiempo }t\\ &Z(t) = \text{ el número de unidades en la zona de alisado en el tiempo }t\\ &S_X = \{0,1,2,3\} \\ &S_Y = \{0,1,2\} \\ &S_Z = \{0,1,2\}&\\ &W(t) = \{(X(t), Y(t), Z(t)\} \\ &S_W = \{(i,j,k), \forall i \in S_X, j \in S_Y, k \in S_Z\} & \end{eqnarray*}\]Vamos a crear el espacio de estados del problema en R, haciendo tres recorridos, de la siguiente forma:
estados = c()
for(i in 0:3){
for(j in 0:2){
for(k in 0:2){
estados = c(estados,paste(i,j,k, sep = ","))
}
}
}
estados## [1] "0,0,0" "0,0,1" "0,0,2" "0,1,0" "0,1,1" "0,1,2" "0,2,0" "0,2,1"
## [9] "0,2,2" "1,0,0" "1,0,1" "1,0,2" "1,1,0" "1,1,1" "1,1,2" "1,2,0"
## [17] "1,2,1" "1,2,2" "2,0,0" "2,0,1" "2,0,2" "2,1,0" "2,1,1" "2,1,2"
## [25] "2,2,0" "2,2,1" "2,2,2" "3,0,0" "3,0,1" "3,0,2" "3,1,0" "3,1,1"
## [33] "3,1,2" "3,2,0" "3,2,1" "3,2,2"length(estados)## [1] 36Podemos ver que nuestra cadena de Markov tiene 36 estados. Ahora bien, tenemos la siguiente formulación general para las tasas de transición entre estados:
\[\begin{equation} r_{i,j}= \begin{cases} 0.05 & \text{si}\ i<3, l=i+1,m=j,n=k \\ (11/12)\cdot 0.2 & \text{si}\ k>0,l=i,m=j,n=k-1 \\ 0.125\cdot min(j,2) & \text{si} \ j>0,k<2,l=i,m=j-1,n=k+1 \\ 0.104 & \text{si} \ i>0,j<2,l=i-1,m=j+1,n=k\\ 0 & \textit{dlc} \end{cases} \end{equation}\]Con la formulación general, podemos definir la matriz de tasas de transición. Primero, vamos a crear una matriz de dimensión \(36 \times 36\) y, para inicializarla, todas las entradas de la misma serán iguales a 0. También vamos a definir los nombres de las filas y columnas de la matriz.
matrizQ=matrix(0, nrow=length(estados), ncol=length(estados))
dimnames(matrizQ) = list(estados, estados)Ahora, utilizando la formulación general, vamos a determinar los valores de cada entrada de la matriz. Nuestro espacio de estados es un vector que contiene elementos de tipo String, y cada elemento contiene información sobre el estado de las variables X(t), Y(t) y Z(t). Por ejemplo, el estado 1 es representado de la siguiente forma:
estados[1]## [1] "0,0,0"Si queremos conocer el valor del estado de cada variable (X(t), Y(t) y Z(t)), debemos separar cada elemento del vector por comas. Esto se realiza mediante la función strsplit(x,split) de R. El parámetro x es un vector que contiene elementos de tipo character o string, y deseamos separar cada uno de estos elementos. El parámetro split es la expresión que se utilizará para separar cada elemento. Vamos a realizar un ejemplo para separar el estado 1 (“0,0,0”):
strsplit(estados[1], ",")## [[1]]
## [1] "0" "0" "0"Conociendo el uso de la función strsplit, podemos hacer dos recorridos sobre el espacio de estados para crear la matriz de transición:
for (fila in estados){
for (columna in estados){
#Para conocer los valores de cada variable de estado, vamos a dividir cada estado en sus tres componentes, utilizando strsplit. Para mayor facilidad, convertimos los estados de cada una de las variables en un valor numérico.
i = as.numeric(strsplit(fila, ",")[[1]][1])
j = as.numeric(strsplit(fila, ",")[[1]][2])
k = as.numeric(strsplit(fila, ",")[[1]][3])
#Una vez más, realizamos la división de cada estado en sus tres componentes.
l = as.numeric(strsplit(columna, ",")[[1]][1])
m = as.numeric(strsplit(columna, ",")[[1]][2])
n = as.numeric(strsplit(columna, ",")[[1]][3])
#Definición de la matriz de tasas de transición, de acuerdo con la formulación general mostrada antes.
#Llegadas a la primera zona (prensas presión)
if(i<3 & l==i+1 & m==j & n==k){
matrizQ[fila,columna] = 0.05
}
#Procesamiento máquina satinadora (zona 3)
if(k>0 & l==i & m==j & n==k-1){
matrizQ[fila,columna] = (11/12)*0.2
}
#Procesamiento zona recubrimiento (zona 2)
if(j>0 & k<2 & l==i & m==j-1 & n==k+1){
matrizQ[fila,columna] = 0.125*min(j,2)
}
#Procesamiento prensa presión (zona 1)
if(i>0 & j<2 & l==i-1 & m==j+1 & n==k){
matrizQ[fila,columna] = 0.1036605
}
}
}En los recorridos anteriores únicamente definimos las tasas de transición entre estados diferentes; no obstante, también debemos asignar los valores correspondientes a la diagonal de la matriz. Como es sabido, el valor de la diagonal será el negativo de la suma de toda la fila:
for(fila in estados){
for(columna in estados){
if(fila==columna){
matrizQ[fila,columna] = -rowSums(matrizQ)[fila]
}
}
}Ahora, vamos a crear la cadena de Markov utilizando el paquete markovchain.
library(markovchain)
cadenaContinua <- new("ctmc", states = estados,
byrow = TRUE, generator = matrizQ)En el literal \(\textbf{b}\) se indica que Gerente de Planta trabaja de 8:00am a 3:00pm. Hoy a las 8:00am él observó que no había ningún rollo siendo procesado en toda el área de fabricación, y desea conocer cuántas veces va a observar los_Input Buffers_ de las tres zonas completamente llenos, hasta las 3:00pm. De 8:00am a 3:00pm hay 7 horas, es decir, 25,200 segundos. Así pues, necesitamos saber cuántas transiciones ocurren hasta 25,200 segundos con el fin de calcular la matriz de ocupación para resolver este literal.
Para esto, necesitamos calcular el valor esperado del tiempo en el que ocurren las transiciones. El tiempo \(T_1\) de la primera transición es:
\[E[T_1]= \alpha \cdot \vec{L}^T\]
El tiempo esperado de las transiciones siguientes se halla de la siguiente forma:
\[E[T_2]= E[T_1]+\alpha \cdot \textbf{P} \cdot \vec{L}^T\] \[â¦\]
\[E[T_k ]= E[T_{k-1}]+\alpha \cdot \textbf{P}^{k-1} \cdot \vec{L}^T\]
Donde \(T_k\) es el tiempo de la k-ésima transición de la cadena a partir del estado inicial.
Para hallar estos tiempos, en primer lugar, creamos un vector \(\vec{L}\) de constantes, tal que su i-ésimo elemento es \(1/r_i\).Cabe mencionar que \(r_i\) es la tasa de permanencia en el estado \(i\), es decir, es igual a \(-q_{i,i}\).
\[\vec{L}=\begin{bmatrix} 1/r_{(0,0,0)} \\ 1/r_{(0,0,1)} \\ ...\\ 1/r_{(3,2,2)}\end{bmatrix}\] Vamos a crear el vector \(\vec{L}\) en R:
L=rep(0,36)
for (i in 1:36){
L[i]=-1/matrizQ[i,i]
}
L## [1] 20.000000 4.285714 4.285714 5.714286 2.790698 4.285714 3.333333
## [8] 2.068966 4.285714 6.507853 2.967413 2.967413 3.588596 2.164531
## [15] 2.967413 3.333333 2.068966 4.285714 6.507853 2.967413 2.967413
## [22] 3.588596 2.164531 2.967413 3.333333 2.068966 4.285714 9.646876
## [29] 3.484395 3.484395 4.373296 2.427221 3.484395 4.000000 2.307692
## [36] 5.454545Ahora, debemos definir el vector de distribución de probabilidades iniciales \(\vec{\alpha}\). Se indica que la cadena inicia sin ningún rollo en proceso,luego:
\[\begin{equation} \alpha_{i,j,k}= \begin{cases} 1 & \text{si}\ (i,j,k)=(0,0,0) \\ 0 & \text{dlc} \end{cases} \end{equation}\]De este modo:
\[\vec{\alpha}=\begin{bmatrix} 1 \\ 0 \\ 0 \\ ...\\ 0 \end{bmatrix}\] Vamos a crear el vector \(\vec{\alpha}\) en R:
#Inicializamos el vector
alpha=rep(0,36)
#Asiganmos un nombre a cada estado
names(alpha)=estados
#Asignamos el valor de 1 al estado inicial
alpha["0,0,0"]=1
alpha## 0,0,0 0,0,1 0,0,2 0,1,0 0,1,1 0,1,2 0,2,0 0,2,1 0,2,2 1,0,0 1,0,1 1,0,2
## 1 0 0 0 0 0 0 0 0 0 0 0
## 1,1,0 1,1,1 1,1,2 1,2,0 1,2,1 1,2,2 2,0,0 2,0,1 2,0,2 2,1,0 2,1,1 2,1,2
## 0 0 0 0 0 0 0 0 0 0 0 0
## 2,2,0 2,2,1 2,2,2 3,0,0 3,0,1 3,0,2 3,1,0 3,1,1 3,1,2 3,2,0 3,2,1 3,2,2
## 0 0 0 0 0 0 0 0 0 0 0 0La matriz \(\textbf{P}\) es la matriz de las probabilidades de transición a un paso de la cadena embebida. Por lo tanto, vamos a hallar la matriz de probabilidades de la cadena embebida:
embebida=matrix(0,nrow=36,ncol=36)
for(i in 1:36){
for (j in 1:36){
if(i!=j){
embebida[i,j]=matrizQ[i,j]/(-matrizQ[i,i])
}
}
}Vamos a calcular el valor esperado del tiempo de las transiciones, hasta la transición 5,000, y almacenaremos estos tiempos en el vector VecT. Ya que calculamos los tiempos para 5,000 transiciones, el vector tendrá 5,000 posiciones, y vamos a inicializarlo con valores de 0.
VecT=rep(0,5000)
#Calculamos el valor esperado del tiempo de la primera transición:
VecT[1]=alpha%*%LPara las siguientes transiciones, vamos a realizar un recorrido. Recordemos que para elevar matrices en R debemos utilizar la librería expm:
library(expm)
for(i in 2:5000){
VecT[i]=VecT[i-1]+alpha%*%(embebida%^%(i-1))%*%L
}El número promedio \(n\) de transiciones que ocurren en la cadena continua en 25,200 segundos está en el intervalo: \(4690\leq n \leq4691\)
VecT[4690]## [1] 25197.75VecT[4691]## [1] 25202.6Utilizamos \(n=4,691\) como aproximación, y hallamos la matriz de ocupación \(\textbf{M}^{4,691}\).
#Se inicializa la matriz con la identidad.
matrizOcupacion=diag(36)
#Nombres de los estados:
dimnames(matrizOcupacion)=list(estados,estados)
#Creación de la matriz
for (i in 1:4691){
matrizOcupacion=matrizOcupacion+embebida%^%i
}El elemento que se solicita en este literal es (0,0,0),(3,2,2) ya que se inicia sin rollos en el área de fabricación, y se desea conocer cuántas veces todas las áreas se encuentran llenas. Obtenemos:
matrizOcupacion["0,0,0","3,2,2"]## [1] 0.6006451Luego, el número de veces es aproximadamente 0.6.
En el literal \(\textbf{c}\), debemos calcular el tiempo promedio que un rollo de papel permanece en el área de fabricación. El tiempo que un rollo permanece en el área (W) se puede hallar mediante la Ley de Little:
\[L=\lambda \cdot W \] Podemos hallar el número de rollos promedio en estado estable (L), así:
\[L=\sum_{(i,j,k) \in S_W} (i+j+k)\cdot \pi_{i,j,k} \] Para esto debemos verificar que la cadena sea irreducible y luego, hallar las probabilidades en estado estable:
is.CTMCirreducible(cadenaContinua)## [1] TRUEVamos a realizar el cálculo de \(L\) de acuerdo con la ecuación presentada anteriormente:
#Inicializamos el valor de L en 0.
LEstadoEstable=0
#Vamos a crear la variable indice: el número del estado que estamos evaluando.
indice=0
for (e in estados){
indice=indice+1
i = as.numeric(strsplit(e, ",")[[1]][1])
j = as.numeric(strsplit(e, ",")[[1]][2])
k = as.numeric(strsplit(e, ",")[[1]][3])
LEstadoEstable=LEstadoEstable+(i+j+k)*steadyStates(cadenaContinua)[indice]
}
LEstadoEstable## [1] 1.457673+0iCon esto, obtenemos que: L=1.458 rollos.
En el literal \(\textbf{d}\), se indica que en este momento hay un rollo en cada una de las zonas del área de fabricación, y se desea conocer cuál es el tiempo promedio que transcurrirá antes de que las tres zonas se encuentren totalmente llenas. Para responder a esta pregunta, debemos hallar el tiempo de primera pasada \(m_{(1,1,1),(3,2,2)}\).
Se hallará este tiempo volviendo el estado (3,2,2) absorbente, y calculando el tiempo antes de absorción dado el estado inicial (1,1,1). Definimos entonces una nueva matriz de tasas de transición, en donde todos los elementos de la fila correspondiente al estado (3,2,2) toman el valor de 0:
matrizAbsorbente=matrizQ
matrizAbsorbente["3,2,2",]=0Para hallar el tiempo antes de absorción debemos calcular \(-\textbf{U}^{-1}\), donde \(\textbf{U}\) es la submatriz de las tasas de transición entre estados transitorios. Creamos entonces la matriz \(\textbf{U}\) así:
estadosTransitorios=estados[-36]
matrizU=matrizAbsorbente[-36,-36]
dimnames(matrizU)=list(estadosTransitorios, estadosTransitorios)Finalmente, calculamos \(-U^(-1)\) y sumamos todos los elementos de la fila correspondiente al estado (1,1,1), dado que este es el estado inicial:
matrizTiemposAbsorcion=-solve(matrizU)
dimnames(matrizTiemposAbsorcion)=list(estadosTransitorios,estadosTransitorios)
sum(matrizTiemposAbsorcion["1,1,1",])## [1] 44362.98Entonces el tiempo de primera pasada al estado (3,2,2) es de 44,362.98 segundos. También es posible realizar este cálculo con la función ExpectedTime(C,i,j) de la librería markovchain. El parámetro C es la cadena de Markov de tipo CTMC, i es el índice del estado incial, y j es el índice del estado final. En este caso, tenemos:
ExpectedTime(cadenaContinua,14,36)## [1] 44362.98