Complementaria 4: Cadenas de Markov en Python III

Cadena de Markov de Tiempo Continuo

En esta complementaria nos concentraremos en la implementación de cadenas de markov en tiempo continuo. Para esto, resolveremos el ejercicio propuesto en el archivo Complementaria 4.pdf

Modele el sistema como una cadena de Markov

En este caso se tienen dos variables. Por un lado, nos interesa guardar el $X(t)$: el nivel de corrosión de la bomba; y $Y(t)$: el nivel de impacto del ambiente.

El espacio de estados de $X(t)$ es $S_X=\{1,\dots,5\}$, el de $Y(t)$ es $S_Y=\{\text{alto},\text{bajo}\}$, y la cadena se puede describir con la variable compuesta $Z(t)=(X(t),Y(t))$ con espacio de estados $S_Z=S_X\times S_Y$.

\[\begin{split} Q_{(i,j)\to(i',j')}= \begin{cases} \lambda_{\text{Alto}}^{-1}, & i'=i+1,\ j'=j,\ i<5,\\ \lambda_{\text{Bajo}}^{-1}, & i'=i-1,\ j'=j,\ i>1,\\ \lambda_{\text{Alto}\to\text{Bajo}}^{-1}, & i'=i,\ j=\text{Alto},\ j'=\text{Bajo},\\ \lambda_{\text{Bajo}\to\text{Alto}}^{-1}, & i'=i,\ j=\text{Bajo},\ j'=\text{Alto},\\ 0, & \text{d.l.c.} \end{cases} \end{split}\]

Con esto, podemos implementar.

#Importar las librerías
import numpy as np
from jmarkov.ctmc import ctmc

#Tasas a utilizar
tasa_ambiente_alto = 0.1
tasa_ambiente_bajo = 0.2
tasa_corrosion_alto_bajo = 0.3
tasa_corrosion_bajo_alto = 0.4

#Crear espacios de estados
estados_corrosion = [i for i in range(1,6)]
estados_ambiente = ['Alto', 'Bajo']

estados = []
for i in estados_corrosion:
    for j in estados_ambiente:
        estado = ",".join((str(i),str(j)))
        estados.append(estado)

estados = np.array(estados)
estados

array(['1,Alto', '1,Bajo', '2,Alto', '2,Bajo', '3,Alto', '3,Bajo',
       '4,Alto', '4,Bajo', '5,Alto', '5,Bajo'], dtype='<U6')

Se define la matriz en ceros y se recorre por fila y por columna para ir ingresando las tasas según las condiciones de la formulación. Dado que estamos en tiempo continuo se tiene que llenar la diagonal para que cada fila sume cero.

#Crear una matriz de ceros

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

#Llenar la matriz de transición

for fila, estado_inicial in enumerate(estados):
    i = int(estado_inicial.split(",")[0])
    j = estado_inicial.split(",")[1]
    for columna, estado_futuro in enumerate(estados):
        i_ = int(estado_futuro.split(",")[0])
        j_ = estado_futuro.split(",")[1]
        if i_==i+1 and j_==j and i<5:
            matrizQ[fila,columna] = tasa_ambiente_alto
        elif i_==i-1 and j_==j and i>1:
            matrizQ[fila,columna] = tasa_ambiente_bajo
        elif i_==i and j=="Alto" and j_=="Bajo":
            matrizQ[fila,columna] = tasa_corrosion_alto_bajo
        elif i_ ==i and j=="Bajo" and j_=="Alto":
            matrizQ[fila,columna] = tasa_corrosion_bajo_alto

# Llenamos la diagonal
for fila, estado_inicial in enumerate(estados):
    matrizQ[fila,fila] = -sum(matrizQ[fila,])

np.sum(matrizQ, axis=1)

#Creamos el objeto cadena

cadena = ctmc(matrizQ)

La segunda pregunta nos dice: Suponga que la bomba inicia operación con el nivel mínimo de corrosión y solo se sabe que, con 50% de probabilidad, el ambiente será de alto impacto. Determine la probabilidad de que, pasada 1 semana de operación, la bomba esté en el nivel máximo de corrosión. Esto es una probabilidad en el transitorio. Como no conocemos con exactitud desde qué estado se inicia se define el vector de probabilidades iniciales de la siguiente manera:

\[\begin{split} \alpha(i,j)=\begin{cases} 0.5, & i=1,\ j\in\{\text{alto},\text{bajo}\},\\ 0, & \text{d.l.c.} \end{cases} \end{split}\]

La probabilidad de que tras 1 semana ($24\times 7$ horas) la bomba esté en el nivel 5 es: $$ \sum_{j\in S_Y} \big[\ \alpha\, e^{Q\,(24\times 7)}\ \big]_{(5,j)}. $$

Para calcular en Python la distribución de probabilidad mencionada se utilizará el método transient_probabilities de la libreria jmarkov que recibe como parametro el vector inicial alpha y el número de pasos necesarios $24\times 7$

#Definir la distribución inicial
alpha = np.zeros(len(estados))
alpha[np.where(estados== "1,Alto")]=0.5
alpha[np.where(estados== "1,Bajo")]=0.5

#Analizar las probabilidades de estado transitorio
probs_transitorio = cadena.transient_probabilities(24*7, alpha)
estado_interes1 = np.where(estados== "5,Alto")
estado_interes2 = np.where(estados== "5,Bajo")

#Calcular la probabilidad de interés
prob_interes = probs_transitorio[estado_interes1]+probs_transitorio[estado_interes2]
prob_interes

array([0.03225762])

Luego, nos preguntan considere el caso en que la máquina inicia su operación en un ambiente con alto impacto y con un nivel de corrosión de 3 (medio). Se observa la máquina a lo largo de 2 semanas. Determine el tiempo esperado (en horas) que pasa la máquina con un nivel de corrosión mayor al estado inicial. Provea una expresión que permita calcular esta cantidad e indique cómo calcularía cada uno de los elementos empleados en su expresión.

Esto se refiere a un tiempo de ocupación para $24\times 7\times 24$ horas, que calculamos como

\[ \sum_{i=4}^{5}\ \sum_{j\in S_Y} \left[ \int_{0}^{2\times7\times24} e^{Q u}\,du \right]_{(3,\text{alto}),(i,j)}. \]

Para calcular en Python el tiempo mencionado se utilizará el método occupation_time de la libreria jmarkov que recibe como parámetro el tiempo de interés $2\times 7\times 24$

#Calcular el tiempo de ocupación
tiempo_ocupacion = cadena.occupation_time(2*7*24)
tiempo_interes = 0
estado_inicial = np.where(estados == "3,Alto")

# Calcular el tiempo de ocupación para estados de interés
for columna, estado_futuro in enumerate(estados):
    i_ = int(estado_futuro.split(",")[0])
    j_ = estado_futuro.split(",")[1]
    if i_>=4:
        tiempo_interes += tiempo_ocupacion[estado_inicial, columna][0][0]

tiempo_interes

34.4620179949542

Finalmente, en el último literal nos piden considerar el caso en que la máquina inicia su operación en un ambiente con bajo impacto y con un nivel de corrosión de 2 (bajo). ¿Cuánto tiempo (en horas) se espera que pase antes de que la máquina llegue a un estado con el máximo nivel de corrosión en un ambiente de alto impacto?

Esto se trata de un tiempo de primera pasada.

La máquina inicia en $Z(0)=(2,\text{bajo})$. El tiempo promedio de primera pasada al estado objetivo $(5,\text{alto})$ se denota $$ m_{(2,\text{bajo})\to(5,\text{alto})}. $$

Sea $\bar{Q}$ la matriz $Q$ sin la fila/columna del estado objetivo. Entonces

\[ m = \bar Q^{-1}\,\mathbf{1},\quad m_{(2,\text{bajo})\to(5,\text{alto})} = [m]_{(2,\text{bajo})}. \]

Equivalente: resolver el sistema lineal

\[ r(i,j)\,m_{(i,j)\to(5,\text{alto})} = 1 + \sum_{(k,\ell)\neq(5,\text{alto})} q_{(i,j),(k,\ell)}\, m_{(k,\ell)\to(5,\text{alto})}. \]

Para calcular en Python el tiempo mencionado se utilizará el método first_time_passage de la libreria jmarkov que recibe como parametro el estado futuro. Además se debe indicar el estado inicial.

# Calcular el tiempo de primera pasada 
#estados de interés
estado_inicial = np.where(estados=='2,Bajo')
estado_futuro = np.where(estados == '5,Alto'
)
# Calcular el tiempo
cadena.first_passage_time(estado_futuro)[estado_inicial][0][0]

276.4706507868981

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