Tarea unidad 2 - Simulación de Procesos y Sistemas - CMTD

Al comienzo de cada semana, el estado de una máquina se determina midiendo la cantidad de corriente eléctrica que utiliza. En función de su lectura de amperaje, la máquina se clasifica en uno de los cuatro estados siguientes: bajo, medio, alto o fallido. Una máquina en estado bajo tiene una probabilidad de \(a\), \(b\) y \(c\) de estar en el estado medio, alto o fallido, respectivamente, al comienzo de la siguiente semana. Una máquina en estado medio tiene una probabilidad de \(d\) y \(e\) de estar en estado alta o fallida, respectivamente, al inicio de la siguiente semana; no puede, por sí sola, pasar al estado bajo. Una máquina en estado alto tiene una probabilidad de \(f\) de estar en el estado fallido al comienzo de la siguiente semana; no puede, por sí misma, pasar al estado bajo o medio. Si una máquina se encuentra en estado de fallo al comienzo de la semana, se inicia inmediatamente la reparación de la máquina para que (con probabilidad 1) esté en el estado bajo al comienzo de la semana siguiente. La máquina funciona 24x7.

Los datos que tienes que utilizar para resolver las preguntas son: probabilidades: (a=0.1;b=0.06;c=0.04;d=0.19;e=0.06;f=0.2) n=2

1. Modeliza este proceso como una CMTD y obtén la matriz de transición de 1 paso. Responde con la suma de los elementos de la diagonal de dicha matriz.

# Hemos de realizar una matriz 4x4 para las probabilidades de transición a n
# pasos, según el enunciado la probabilidad a de pasar de bajo a medio, b de
# pasar a alto y c de pasar a fallido, la probabilidad de pasar de medio a alto
# es de d, y a fallido es de e, la prob. de pasar de alto a fallido es de f.
# Por último, si la máquina falla en el siguiente paso comienza automáticamente
# en bajo. Por lo cual la matriz de transición a 1 paso quedaría así: 
library(markovchain)
## Package:  markovchain
## Version:  0.9.1
## Date:     2023-01-20
## BugReport: https://github.com/spedygiorgio/markovchain/issues
estados <- c("bajo", "medio", "alto", "fallido")
a <- 0.1; b <- 0.06; c <- 0.04; d <- 0.19; e <- 0.06; f <- 0.2;
matriz <- matrix(c(1-a-b-c, a, b, c,0, 1-d-e, d, e, 0, 0, 1-f, f, 1, 0, 0, 0),
                 nrow = 4, ncol = 4, byrow = TRUE,
                 dimnames = list(estados, estados))

# La prob. de pasar de bajo a bajo es 1 - las probabilidades acumuladas de
# pasar a medio, alto y fallo respectivamente(0.1 + 0-06 + 0.04). La prob. de
# pasar de medio a medio es de 1 - las probabilidades acumuladas de pasar a alto
# y fallo respectivamente (0.19 + 0.06). La prob. de pasar de alto a alto es de
# 1 - la probabilidad de pasar a fallo (0.2), y si se encuentra en estado de
# fallo la semana siguiente pasará automáticamente al estado bajo (p = 1).

suma_diagonal <- sum(diag(matriz))
print(matriz)
##         bajo medio alto fallido
## bajo     0.8  0.10 0.06    0.04
## medio    0.0  0.75 0.19    0.06
## alto     0.0  0.00 0.80    0.20
## fallido  1.0  0.00 0.00    0.00
cat("\nLa suma de los elementos de la diagonal es:", suma_diagonal)
## 
## La suma de los elementos de la diagonal es: 2.35
proceso <- new("markovchain", states = estados, 
               byrow = TRUE, transitionMatrix = matriz)
proceso
## Unnamed Markov chain 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  bajo, medio, alto, fallido 
##  The transition matrix  (by rows)  is defined as follows: 
##         bajo medio alto fallido
## bajo     0.8  0.10 0.06    0.04
## medio    0.0  0.75 0.19    0.06
## alto     0.0  0.00 0.80    0.20
## fallido  1.0  0.00 0.00    0.00
plot(proceso)

2. Si una máquina nueva siempre comienza en estado bajo, ¿cuál es la probabilidad de que la máquina esté en estado de fallo después de 3 semanas?

# Definimos la función para generar la matriz de transición.
# Matriz de probs. transición de n pasos
ptran.n=function(ptran,n){
  # ptran es la matriz de transición de 1 paso
  # n son los pasos a dar
  i=1
  p=ptran
  while(i<n){
    p=p%*%ptran
    i=i+1
  }
return(p)  
}
# Debemos multiplicar el vector del estado inicial por la matriz de transición correspondiente.
# En este caso comenzamos en estado bajo, y el estado final es fallo después de 3 semanas.
est_inicial_2 <- c(1,0,0,0)
# matriz de transición de 3 pasos
matriz_3 <- ptran.n(matriz,3); matriz_3
##          bajo    medio     alto fallido
## bajo    0.594 0.184250 0.162250 0.05950
## medio   0.131 0.427875 0.346075 0.09505
## alto    0.320 0.020000 0.524000 0.13600
## fallido 0.680 0.155000 0.115000 0.05000
cat("\n La probabilidad de pasar del estado bajo a fallo en 3 semanas es de", matriz_3[1,4])
## 
##  La probabilidad de pasar del estado bajo a fallo en 3 semanas es de 0.0595
cat("\n Aproximadamente un 5.95%.")
## 
##  Aproximadamente un 5.95%.

3. Y si hoy arranca en estado bajo, ¿cuántas horas (en promedio) transcurrirán hasta que dé un fallo?

meanFirstPassageTime(proceso, "fallido")
##  bajo medio  alto 
##  10.4   7.8   5.0
cat("\n En promedio, si comienza en estado bajo transcurrirán 10.4 transiciones hasta que de un fallo por primera vez.")
## 
##  En promedio, si comienza en estado bajo transcurrirán 10.4 transiciones hasta que de un fallo por primera vez.
cat("\n Como las transiciones son semanales, para conseguir el promedio en horas multiplicaremos el número de transiciones por 168.")
## 
##  Como las transiciones son semanales, para conseguir el promedio en horas multiplicaremos el número de transiciones por 168.
cat("\n Necesitando así ", 10.4*168, "horas transcurridas hasta que de un fallo.")
## 
##  Necesitando así  1747.2 horas transcurridas hasta que de un fallo.
cat("\n Si comenzamos en estado medio, necesitaremos ", 7.8*168, "horas y si comenzamos en alto necesitaremos en promedio ", 5*168, "horas.")
## 
##  Si comenzamos en estado medio, necesitaremos  1310.4 horas y si comenzamos en alto necesitaremos en promedio  840 horas.

4. Si hoy está en modo fallo, ¿cuántas horas (en promedio) transcurrirán hasta que vuelva a dar un fallo?

mFPT_FF <- meanFirstPassageTime(proceso);mFPT_FF
##         bajo medio      alto fallido
## bajo     0.0    14 10.764706    10.4
## medio    8.8     0  6.823529     7.8
## alto     6.0    20  0.000000     5.0
## fallido  1.0    15 11.764706     0.0
mFPT_FF[4,4]
## [1] 0
cat("\n Si comenzamos en estado fallido, ya estamos en estado de fallo, por lo que han de transcurrir 0 horas,")
## 
##  Si comenzamos en estado fallido, ya estamos en estado de fallo, por lo que han de transcurrir 0 horas,

Si la pregunta se refiere al tiempo te recurrencia y no al de primer paso, lo podemos calcular así:

MRT <- meanRecurrenceTime(proceso); MRT
##      bajo     medio      alto   fallido 
##  2.280000  5.700000  3.352941 11.400000
cat("\n El tiempo te recurrencia (en promedio) es de", MRT[4],"transiciones, unas", MRT[4]*168,"horas.")
## 
##  El tiempo te recurrencia (en promedio) es de 11.4 transiciones, unas 1915.2 horas.

5. Si durante esta semana la máquina trabaja en estado bajo, transcurridas \(n\) semanas más, ¿cuántos días habrá estado funcionando la máquina?

(Ayuda: Verifica que el número total de semanas en funcionamiento va ser de \(n+1\)).

# Dado que n=2, debemos comprobar cuanto tiempo ha estado en funcionamiento la máquina durante 3 semanas (n+1).
set.seed(123)
state_seq <- markovchainSequence(2, markovchain = proceso, t0="bajo", include.t0 = TRUE);state_seq
## [1] "bajo" "bajo" "bajo"
cat("\n Según la simulación, comenzando en estado bajo las dos semanas siguientes continuará en estado bajo.")
## 
##  Según la simulación, comenzando en estado bajo las dos semanas siguientes continuará en estado bajo.
cat("\n Lo que hace un total de", sum(state_seq == "bajo", state_seq == "medio", state_seq == "alto"), "semanas, o", sum(state_seq == "bajo", state_seq == "medio", state_seq == "alto")*7, "días.")
## 
##  Lo que hace un total de 3 semanas, o 21 días.

6. A lo largo de un año, ¿cuánto tiempo, en porcentaje, estará funcionando la máquina?

set.seed(123)
state_seq_anual <- markovchainSequence(52, markovchain = proceso);
semanas_trabajadas <- sum(state_seq_anual == "bajo", state_seq_anual == "medio", state_seq_anual == "alto")
cat("\n La máquina estará funcionando un ", (semanas_trabajadas*100)/52, "% del año.")
## 
##  La máquina estará funcionando un  90.38462 % del año.

7. A lo largo de un año, ¿cuántas semanas estará operativa la máquina? (52 semanas en 1 año)

# Utilizando el mismo código (ya que sigue siendo una simulación a 1 año):
set.seed(123)
state_seq_anual <- markovchainSequence(52, markovchain = proceso);
semanas_trabajadas <- sum(state_seq_anual == "bajo", state_seq_anual == "medio", state_seq_anual == "alto")
cat("\n La máquina estará funcionando", semanas_trabajadas, " semanas durante el año.")
## 
##  La máquina estará funcionando 47  semanas durante el año.

8. ¿Cuál es el beneficio semanal a largo plazo?

Cada semana que la máquina está en estado bajo, se obtiene un beneficio de 1000 euros; cada semana que la máquina está en el estado medio, se obtiene un beneficio de 500 euros; cada semana que la máquina está en estado alto, se obtiene un beneficio de 400 euros; y la semana en la que hay un fallo, se incurre en un coste de 700 euros por la reparación.

# Utilizaremos un periodo de 10 años (520 semanas) para la simulación.
set.seed(123)
state_seq_trianual <- markovchainSequence(520, markovchain = proceso);
semanas_bajo <- sum(state_seq_trianual == "bajo")
semanas_medio <- sum(state_seq_trianual == "medio")
semanas_alto <- sum(state_seq_trianual == "alto")
semanas_fallo <- sum(state_seq_trianual == "fallido")
semanas_trabajadas_trianual <- sum(semanas_bajo + semanas_medio + semanas_alto)
ingresos_trianual <- 1000*semanas_bajo + 500*semanas_medio + 400*semanas_alto
gastos_trianual <- 700*semanas_fallo
total_trianual <- ingresos_trianual - gastos_trianual
cat("\n La máquina estará funcionando", semanas_trabajadas_trianual, " semanas durante el año.")
## 
##  La máquina estará funcionando 473  semanas durante el año.
cat("\n Estará funcionando", semanas_bajo, "semanas en estado bajo,", semanas_medio, "semanas en estado medio", semanas_alto, "semanas en estado alto y fallará durante", semanas_fallo, "semanas.")
## 
##  Estará funcionando 257 semanas en estado bajo, 68 semanas en estado medio 148 semanas en estado alto y fallará durante 47 semanas.
cat("\n Con unos beneficios de 1000 euros por cada semana funcionando en estado bajo, de 500 euros en medio, de 400 euros en estado alto, y un coste de reparación de 700 euros cada vez que falla:")
## 
##  Con unos beneficios de 1000 euros por cada semana funcionando en estado bajo, de 500 euros en medio, de 400 euros en estado alto, y un coste de reparación de 700 euros cada vez que falla:
cat("\n Obtendremos unos beneficios de:", ingresos_trianual, "euros, y unos gastos de", gastos_trianual, ". Haciendo así un total de", total_trianual, "euros de beneficios durante los 3 años.")
## 
##  Obtendremos unos beneficios de: 350200 euros, y unos gastos de 32900 . Haciendo así un total de 317300 euros de beneficios durante los 3 años.
cat("\n Siendo esto un beneficio semanal de", total_trianual/520, " euros semanales.")
## 
##  Siendo esto un beneficio semanal de 610.1923  euros semanales.

9. ¿Merece la pena esta nueva política de mantenimiento en términos de beneficios esperados (sí/no)?

Se ha sugerido cambiar la política de mantenimiento de la máquina. Si al comienzo de una semana la máquina está en estado alto, la máquina se deja fuera de servicio y es reparada para que al inicio de la siguiente semana vuelva a estar en el estado bajo. Cuando se realiza una reparación se incurre en un coste de 600 euros. Registra además, cuál es el beneficio semanal que se conseguiría.

# Al entrar en estado alto se pone fuera de servicio y la semana siguiente tomará
# el estado bajo, por lo que adopta las probabilidades que tenía en la matriz
# anterior el estado fallo, esto es: 

matriz_v2 <- matrix(c(1-a-b-c, a, b, c,0, 1-d-e, d, e, 1, 0, 0, 0, 1, 0, 0, 0),
                 nrow = 4, ncol = 4, byrow = TRUE,
                 dimnames = list(estados, estados))

proceso2 <- new("markovchain", states = estados, 
               byrow = TRUE, transitionMatrix = matriz_v2)
proceso2
## Unnamed Markov chain 
##  A  4 - dimensional discrete Markov Chain defined by the following states: 
##  bajo, medio, alto, fallido 
##  The transition matrix  (by rows)  is defined as follows: 
##         bajo medio alto fallido
## bajo     0.8  0.10 0.06    0.04
## medio    0.0  0.75 0.19    0.06
## alto     1.0  0.00 0.00    0.00
## fallido  1.0  0.00 0.00    0.00
plot(proceso2)

# Una vez creada la nueva matriz de transición, y la nueva cadena de Markov,
# procedemos a realizar la simulación y cálculo de los costes actualizados a la
# nueva política de mantenimientos.
# Volvemos a utilizar un periodo de 10 años para comparar los resultados con los
# del ejercicio anterior.
set.seed(123)
state_seq_trianual2 <- markovchainSequence(520, markovchain = proceso2);
semanas_bajo2 <- sum(state_seq_trianual2 == "bajo")
semanas_medio2 <- sum(state_seq_trianual2 == "medio")
semanas_alto2 <- sum(state_seq_trianual2 == "alto")
semanas_fallo2 <- sum(state_seq_trianual2 == "fallido")
semanas_trabajadas_trianual2 <- sum(semanas_bajo2 + semanas_medio2)
ingresos_trianual2 <- 1000*semanas_bajo2 + 500*semanas_medio2
gastos_trianual2 <-  + 600*semanas_alto2 + 600*semanas_fallo2
total_trianual2 <- ingresos_trianual2 - gastos_trianual2
cat("\n La máquina estará funcionando", semanas_trabajadas_trianual2, " semanas durante el año.")
## 
##  La máquina estará funcionando 458  semanas durante el año.
cat("\n Estará funcionando", semanas_bajo2, "semanas en estado bajo,", semanas_medio2, "semanas en estado medio.")
## 
##  Estará funcionando 345 semanas en estado bajo, 113 semanas en estado medio.
cat("\n Se alcanzará el estado alto", semanas_alto2, "semanas con lo cual tendremos que ponerla fuera de servicio, y fallará durante", semanas_fallo2, "semanas.")
## 
##  Se alcanzará el estado alto 34 semanas con lo cual tendremos que ponerla fuera de servicio, y fallará durante 28 semanas.
cat("\n Con unos beneficios de 1000 euros por cada semana funcionando en estado bajo, de 500 euros en medioy unos costes de reparación de 600 euros cada vez que alcanza el estado alto o falla:")
## 
##  Con unos beneficios de 1000 euros por cada semana funcionando en estado bajo, de 500 euros en medioy unos costes de reparación de 600 euros cada vez que alcanza el estado alto o falla:
cat("\n Obtendremos unos beneficios de:", ingresos_trianual2, "euros, y unos gastos de", gastos_trianual2, ". Haciendo así un total de", total_trianual2, "euros de beneficios durante los 3 años.")
## 
##  Obtendremos unos beneficios de: 401500 euros, y unos gastos de 37200 . Haciendo así un total de 364300 euros de beneficios durante los 3 años.
cat("\n Siendo esto un beneficio semanal de", total_trianual2/520, " euros semanales.")
## 
##  Siendo esto un beneficio semanal de 700.5769  euros semanales.
cat("\n El cambio en la política de mantenimiento supone un aumento en los beneficios semanales de", (total_trianual2/520)-(total_trianual/520), "euros, por lo cual merece la pena la nueva política.")
## 
##  El cambio en la política de mantenimiento supone un aumento en los beneficios semanales de 90.38462 euros, por lo cual merece la pena la nueva política.