Complementaria 6: Análisis de tiempos promedio

Complementaria 6: Análisis de tiempos promedio

#

En este tutorial se abordará el tema de análisis de tiempos en cadenas de Markov a través de un ejemplo. Para esto, utilizaremos la librería jmarkov.

Modelación del proceso de producción en una fábrica de papel

Vamos a resolver el problema del archivo “Complementaria 6(Q).pdf” que se encuentra en Bloque Neón.

Literal a

En primer lugar, se modelará la situación como una cadena de Markov de tiempo continuo. Se definen tres variables de estado:

\[ \begin{align}\begin{aligned}\begin{split} X(t) = \text{número de unidades en la zona de prensas en el tiempo } t\\\end{split}\\\begin{split}Y(t) = \text{número de unidades en la zona de recubrimiento en el tiempo } t\\\end{split}\\Z(t) = \text{número de unidades en la zona de alisado en el tiempo } t \end{aligned}\end{align} \]
\[ \begin{align}\begin{aligned}\begin{split} S_X = \{0,1,2,3\}\\\end{split}\\\begin{split}S_Y = \{0,1,2\}\\\end{split}\\\begin{split}S_Z = \{0,1,2\}\\\end{split}\\\begin{split}W(t) = (X(t),Y(t),Z(t))\\\end{split}\\S_W = \{(i,j,k), \forall i \in S_X , j \in S_Y , k \in S_Z\} \end{aligned}\end{align} \]

Vamos a crear el espacio de estados del problema en Python, haciendo tres recorridos, de la siguiente forma:

estados = []

for i in range(4):
    for j in range(3):
        for k in range(3):
            estados.append(f"{i},{j},{k}")
            
estados

['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,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',
 '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']

len(estados)

36

Podemos 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{split} \mathbb{Q}_{(i,j,k) \rightarrow (l,m,n)} = \begin{cases} 0.05 & l=i+1 & m=j & n=k & i<3 \\ 11/12 \cdot 0.2 & l=i & m=j & n=k-1 & k>0 \\ 0.125 \cdot \min(j,2) & l=i & m=j-1 & n=k+1 & j>0 & k<2 \\ 0.104 & l=i-1 & m=j+1 & n=k & i>0 & j<2 \\ 0 & \text{d.l.c} \end{cases} \end{split}\]

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.

import numpy as np

matrizQ = np.zeros((len(estados), len(estados)), dtype = float)

Ahora, vamos a convertir esta matriz en un Dataframe para asignarle nombre a cada fila y columna de nuestra matriz.

import pandas as pd

matrizQ = pd.DataFrame(matrizQ, index=estados, columns=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[0]

'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 texto.split() de Python. El objeto texto es una cadena de caracteres que contiene los elementos de tipo character o string, y deseamos separar cada uno de estos elementos. Dentro de la función splitdebemos indicar cuál será el caracter encargado de separar los elementos (en este caso es la coma). Vamos a realizar un ejemplo para separar el estado 1 (“0,0,0”):

np.array(estados[0].split(','), dtype=int)

array([0, 0, 0])

Conociendo el uso de la función split, podemos hacer dos recorridos sobre el espacio de estados para crear la matriz de transición.

for fila in estados:
    for columna in estados:
        i,j,k  = np.array(fila.split(','), dtype=int)
        
        l,m,n = np.array(columna.split(','), dtype=int)
        
        #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 and l==i+1 and m==j and n==k:
            matrizQ.loc[fila,columna] = 0.05
        
        #Procesamiento máquina satinadora (zona 3)
        if k>0 and l==i and m==j and n==k-1:
            matrizQ.loc[fila,columna] = (11/12)*0.2
        
        #Procesamiento zona recubrimiento (zona 2)
        if j>0 and k<2 and l==i and m==j-1 and n==k+1:
            matrizQ.loc[fila,columna] = 0.125*min(j,2)
            
        #Procesamiento prensa presión (zona 1)
        if i>0 and j<2 and l==i-1 and m==j+1 and n==k:
            matrizQ.loc[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 i in range(len(estados)):
    # Diagonal: negativo de la suma de la fila
    matrizQ.iloc[i, i] = -np.sum(matrizQ.iloc[i, :])

Ahora, vamos a crear la cadena de Markov utilizando la librería jmarkov. Para ello, debemos volver a hacer que la matriz que creamos sea de tipo np.array.

matrizQ = np.array(matrizQ)

from jmarkov.ctmc import ctmc

cadenaContinua = ctmc(matrizQ)

Literal b

En el literal b se indica que el Gerente de Planta trabaja de 8:00 am a 3:00 pm. Hoy a las 8:00 am él observó que no había ningún rollo siendo procesado en toda el área de fabricación, y desea conocer qué proporción del tiempo va a observar los Input Buffers de las tres zonas completamente llenos, hasta las 3:00 pm. De 8:00 am a 3:00 pm hay 7 horas, es decir, 25,200 segundos.

Así pues, necesitamos calcular la matriz de ocupación para resolver este literal. Con la matriz de ocupación podemos estimar cuánto tiempo estarán las tres zonas completamente llenas, es decir, cuánto tiempo durante esas 7 horas el sistema estará en el estado \((3,2,2)\).

Para ello, usaremos la función occupation_time(t) de la librería jmarkov, en donde el argumento t representa el tiempo para el cual se quiere hacer el análisis de la matriz de tiempos de ocupación.

t_total = 25200
matrizTiempos = cadenaContinua.occupation_time(t_total)

Una vez generada la matriz de tiempos de ocupación, es necesario recuperar la fila asociada al estado (0,0,0) y la columna asociada al estado (1,1,1). Para ello, se usará la función index(estado) la cual recupera la posición en la matriz asociada al estadode interés.

tiempo = matrizTiempos[estados.index('0,0,0'), estados.index('3,2,2')]
tiempo

3.279451538778562

Luego, el tiempo es de aproximadamente 3.28 segundos. Sin embargo, debido a que nos preguntan sobre la proporción del tiempo que todas las estaciones permanecen ocupadas, se calcula la proporción dividiendo el tiempo obtenido sobre el tiempo total de observación:

proporcion = tiempo/t_total
proporcion

0.00013013696582454612

Finalmente, obtenemos que la proporción del tiempo que las estaciones permanecen llenas es del 0.013%.

Literal c

En el literal 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 (usando la función is_irreducible()) y luego, hallar las probabilidades en estado estable (con la función steady_state()):

cadenaContinua.is_irreducible()

True

Vamos 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

#Creamos el vector de Estado Estable
estadoEstable = cadenaContinua.steady_state()

for e in estados:
    i,j,k = np.array(e.split(','), dtype=int)
    LEstadoEstable = LEstadoEstable + (i+j+k)*estadoEstable[indice]
    
    indice = indice+1
  
LEstadoEstable

1.4576733281181173

Con esto, obtenemos que: \(L=1.458\) rollos.

Literal d

En el literal 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)}\).

Es posible realizar este cálculo con la función first_passage_time(target)[inicial] de la librería jmarkov. Donde target es el estado futuro al que se quiere llegar. Esta función genera un array, en donde cada elemento corresponde al posible estado inicial sobre el cual se quiere hacer el cálculo de tiempo de primera pasada. De esta manera, el parámetro target debe indicar el estado al que se quiere llegar, y sobre el array resultado se debe llamar la posición asociada al estado (1,1,1) (que corresponde al argumento [inicial] de la explicación anterior). Por lo tanto:

cadenaContinua.first_passage_time(estados.index('3,2,2'))[estados.index('1,1,1')]

array([44362.97918331])

El tiempo promedio que transcurre desde el estado (1,1,1) hasta visitar por primera vez el estado (3,2,2) es de 44362 segundos, o 12.32 horas.

Universidad de los Andes | Vigilada Mineducación. Reconocimiento como Universidad: Decreto 1297 del 30 de mayo de 1964. Reconocimiento personería jurídica: Resolución 28 del 23 de febrero de 1949 Minjusticia. Departamento de Ingeniería Industrial Carrera 1 Este No. 19 A 40 Bogotá, Colombia Tel. (57.1) 3324320 | (57.1) 3394949 Ext. 2880 /2881 http://industrial.uniandes.edu.co