Tarea Unidad 04. Cadenas de Markov de tiempo continuo
Problema planteado
Una empresa realiza un proceso de fabricación con 6 máquinas que trabajan diáriamente sin parada durante 8 horas, 5 días a la semana (de lunes a viernes). Cuando el sistema está en pleno funcionamiento, la empresa genera 200€ de beneficio a la hora, sin embargo, cuando se avería 1 máquina, el beneficio se reduce a la mitad, 100€ a la hora y en el caso de que se averíe una segunda máquina pasaría a tan solo 50€ de beneficio la hora. Ahora bien, cuando se averían 3 o más máquinas el proceso se detiene por completo y no hay beneficios.
El tiempo de fallo de cualquier máquina se distribuye según una distribución de Poisson de media 7 días. El tiempo de reparación por otro lado, es una exponencial con tasa mu.
Objetivo de la fábrica
- El objetivo del análisis es obtener unas conclusiones acerca de la contratación que necesita la empresa para poder producir el máximo beneficio posible.
Metodología
Hemos utilizado la Metodología de Cadenas de Markov de tiempo continuo, por sus siglas (CTMC), construyendo un proceso de nacimiento-muerte, donde el nacimiento es el fallo de una de las 6 máquinas y la muerte es la reparación de la misma.
Para esta metodología es clave conocer los valores de \(\lambda\) y \(\mu\), además, hemos utilizado la función definida como tiempos.ocupación y que necesida de inputs tales como Rmat= matriz de tasas, Ts= tiempo de ocupación, y cal= forma de obtener valores (siempre 1 en este caso). Además hemos utilizado la librería markovchain.
Cuestiones planteadas
1. La empresa necesita saber cúantos reparadores tiene que tener y cúanto debe pagarle a cada uno de ellos, para ello, examinaremos el tiempo de fallo de las máquinas. Para ello estudiaremos cúanto tiempo tarda en pararse la producción, es decir, en que fallen las 6 máquinas, si todas estaban operativas.
estados=as.character(0:6)
nestados=length(estados)
lambda=1/(7*24)
set.seed(2)
unif=runif(1,2,24)
mu=0
R = matrix(nrow = nestados, ncol = nestados, data = 0,
dimnames=list(estados,estados))
R[1,2]= 6*lambda
R[2,1]= mu
R[2,3]= 5*lambda
R[3,2]= 2*mu
R[3,4]= 4*lambda
R[4,3]= 3*mu
R[4,5]= 3*lambda
R[5,4]= 4*mu
R[5,6]= 2*lambda
R[6,5]= 5*mu
R[6,7]= lambda
R[7,6]= 6*mu
Q=R-diag(unlist(apply(R,1,sum)))
maquinas=new("ctmc",states=estados,byrow=TRUE, generator=Q, name="maquinas")
et=ExpectedTime(maquinas,1,7)
cat("\n Sin ningún reparador, la maquinaria tardará",et,"horas en pararse por completo")##
## Sin ningún reparador, la maquinaria tardará 411.6 horas en pararse por completo- Hemos visto que la empresa sin reparadores podría estar produciendo tan solo 411.6 horas, por lo que no es nada factible no contratar reparadores.
Para analizar la viabilidad de beneficio/coste al emplear trabajadores vamos a plantear 3 escenarios diferentes:
- Un trabajador
- Dos trabajadores
- Tres trabajadores
estados=as.character(0:6)
nestados=length(estados)
lambda=1/(7*24)
set.seed(2)
unif=runif(1,2,24)
mu=1/unif
# Matriz de costes
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)
rlimit=ifelse(cal == 1, max(ris),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))
}Un reparador
Contratamos al reparador cuando hay una máquina averiada. ¿Cúanto trabajo va a tener el reparador el primer día de trabajo?¿Y la primera semana?
La contratación de los dos reparadores la realizamos a un mes, por lo que queremos saber el beneficio que nos reportan si ambos cobran 1000€ al mes. (Tendremos en cuenta que un mes tiene 4 semanas laborables).
# Generamos la matriz de tasas, la que muestra la transición del punto i al punto j
# Proceso de nacimiento y muerte, en el que los nacimientos implican averías y las muertes reparaciones
R1 = matrix(nrow = nestados, ncol = nestados, data = 0,
dimnames=list(estados,estados))
R1[1,2]= 6*lambda
R1[2,1]= mu
R1[2,3]= 5*lambda
R1[3,2]= 1*mu
R1[3,4]= 4*lambda
R1[4,3]= 1*mu
R1[4,5]= 3*lambda
R1[5,4]= 1*mu
R1[5,6]= 1*lambda
R1[6,5]= 1*mu
R1[6,7]= lambda
R1[7,6]= 1*mu
# Generamos la matriz generadora
Q1=R1-diag(unlist(apply(R1,1,sum)))
maquinas1=new("ctmc",states=estados,byrow=TRUE, generator=Q1, name="maquinas")
## Primer día
g1.1=tiempos.ocupacion(R1,8,1)
trabajo.primer.dia1=8-(g1.1[2,1])
cat("\n Tiempo en el que un reparador está ocupado durante el primer día de trabajo",round(trabajo.primer.dia1,3),"horas,es decir,un",trabajo.primer.dia1/8*100,"% del tiempo")##
## Tiempo en el que un reparador está ocupado durante el primer día de trabajo 4.849 horas,es decir,un 60.6125 % del tiempo## Primera semana
g1.2=tiempos.ocupacion(R1,8*5,1)
trabajo.primera.semana1=8*5-(g1.2[2,1])
cat("\n Tiempo en el que un reparador está ocupado en una semana",round(trabajo.primera.semana1,3),"horas,es decir,un",trabajo.primera.semana1/(8*5)*100,"% del tiempo")##
## Tiempo en el que un reparador está ocupado en una semana 12.653 horas,es decir,un 31.63325 % del tiempo## Tiempo ocupación un mes para calcular beneficios
beneficios=matrix( c(200,100,50,0,0,0,0),ncol=1)
g1.3=tiempos.ocupacion(R1,5*8*4,1)
finmes1=g1.3%*%beneficios
trab1=matrix(c(1000,1000,1000,1000,1000,1000,1000),ncol=1)
benef.finmes1=finmes1-trab1
cat("\n Un mes aún no se considera estado estacionario, por lo que si influye por qué máquina se comience, por ello diremos que la empresa generará a fin de mes contando un sueldo de 1000€, un beneficio entre:",max(benef.finmes1),"y",min(benef.finmes1))##
## Un mes aún no se considera estado estacionario, por lo que si influye por qué máquina se comience, por ello diremos que la empresa generará a fin de mes contando un sueldo de 1000€, un beneficio entre: 27522.62 y 21439.29Dos reparadores
Contratamos a los dos reparadores cuando se han roto dos máquinas. ¿Cómo de ocupados van a estar los dos trabajadores el primer día?¿Y la primera semana?
La contratación de los dos reparadores la realizamos a un mes, por lo que queremos saber el beneficio que nos reportan si ambos cobran 1000€ al mes. (Tendremos en cuenta que un mes tiene 4 semanas laborables).
# Generamos la matriz de tasas, la que muestra la transición del punto i al punto j
# Proceso de nacimiento y muerte, en el que los nacimientos implican averías y las muertes reparaciones
R2 = matrix(nrow = nestados, ncol = nestados, data = 0,
dimnames=list(estados,estados))
R2[1,2]= 6*lambda
R2[2,1]= mu
R2[2,3]= 5*lambda
R2[3,2]= 2*mu
R2[3,4]= 4*lambda
R2[4,3]= 2*mu
R2[4,5]= 3*lambda
R2[5,4]= 2*mu
R2[5,6]= 2*lambda
R2[6,5]= 2*mu
R2[6,7]= lambda
R2[7,6]= 2*mu
# Generamos la matriz generadora
Q2=R2-diag(unlist(apply(R2,1,sum)))
maquinas2=new("ctmc",states=estados,byrow=TRUE, generator=Q2, name="maquinas")
## Primer día
g2.1=tiempos.ocupacion(R2,8,1)
trabajo.primer.dia=8-(g2.1[3,1]+g2.1[3,2])
cat("\n Tiempo de los dos trabajadores reparando al mismo tiempo el primer día",round(trabajo.primer.dia,3),"horas, es decir,un",trabajo.primer.dia/8*100,"% del tiempo")##
## Tiempo de los dos trabajadores reparando al mismo tiempo el primer día 3.103 horas, es decir,un 38.79 % del tiempo## Primera semana
g2.2=tiempos.ocupacion(R2,8*5,1)
trabajo.primera.semana=8*5-(g2.2[3,1]+g2.2[3,2])
cat("\n Tiempo de los dos trabajadores reparando al mismo tiempo en una semana",round(trabajo.primera.semana,3),"horas, es decir,un",trabajo.primera.semana/(8*5)*100,"% del tiempo")##
## Tiempo de los dos trabajadores reparando al mismo tiempo en una semana 4.214 horas, es decir,un 10.53475 % del tiempo## Tiempo ocupación un mes para calcular beneficios
beneficios=matrix( c(200,100,50,0,0,0,0),ncol=1)
g2.3=tiempos.ocupacion(R2,5*8*4,1)
finmes2=g2.3%*%beneficios
trab2=matrix(c(2000,2000,2000,2000,2000,2000,2000),ncol=1)
benef.finmes2=finmes2-trab2
cat("\n Un mes aún no se considera estado estacionario, por lo que si influye por qué máquina se comience, por ello diremos que la empresa generará a fin de mes contando dos sueldos de 1000€ cada uno, un beneficio entre:",max(benef.finmes2),"y",min(benef.finmes2))##
## Un mes aún no se considera estado estacionario, por lo que si influye por qué máquina se comience, por ello diremos que la empresa generará a fin de mes contando dos sueldos de 1000€ cada uno, un beneficio entre: 26897.28 y 23654.89Tres reparadores
Contratamos a los tres reparadores cuando se han roto tres máquinas. ¿Cómo de ocupados van a estar los dos trabajadores el primer día?¿Y la primera semana?
La contratación de los dos reparadores la realizamos a un mes, por lo que queremos saber el beneficio que nos reportan si ellos cobran 1000€ al mes. (Tendremos en cuenta que un mes tiene 4 semanas laborables).
# Generamos la matriz de tasas, la que muestra la transición del punto i al punto j
# Proceso de nacimiento y muerte, en el que los nacimientos implican averías y las muertes reparaciones
R3 = matrix(nrow = nestados, ncol = nestados, data = 0,
dimnames=list(estados,estados))
R3[1,2]= 6*lambda
R3[2,1]= mu
R3[2,3]= 5*lambda
R3[3,2]= 2*mu
R3[3,4]= 4*lambda
R3[4,3]= 3*mu
R3[4,5]= 3*lambda
R3[5,4]= 3*mu
R3[5,6]= 2*lambda
R3[6,5]= 3*mu
R3[6,7]= lambda
R3[7,6]= 3*mu
# Generamos la matriz generadora
Q3=R3-diag(unlist(apply(R3,1,sum)))
maquinas3=new("ctmc",states=estados,byrow=TRUE, generator=Q3, name="maquinas")
## Primer día
g3.1=tiempos.ocupacion(R3,8,1)
trabajo.primer.dia=8-(g3.1[4,1]+g3.1[4,2]+g3.1[4,3])
cat("\n Tiempo de los dos trabajadores reparando al mismo tiempo el primer día",round(trabajo.primer.dia,3),"horas, es decir,un",trabajo.primer.dia/8*100,"% del tiempo")##
## Tiempo de los dos trabajadores reparando al mismo tiempo el primer día 2.138 horas, es decir,un 26.72 % del tiempo## Primera semana
g3.2=tiempos.ocupacion(R3,8*5,1)
trabajo.primera.semana=8*5-(g3.2[4,1]+g3.2[4,2]+g3.2[4,3])
cat("\n Tiempo de los dos trabajadores reparando al mismo tiempo en una semana",round(trabajo.primera.semana,3),"horas, es decir,un",trabajo.primera.semana/(8*5)*100,"% del tiempo")##
## Tiempo de los dos trabajadores reparando al mismo tiempo en una semana 2.29 horas, es decir,un 5.724 % del tiempo## Tiempo ocupación un mes
beneficios=matrix( c(200,100,50,0,0,0,0),ncol=1)
g3.3=tiempos.ocupacion(R3,5*8*4,1)
finmes3=g3.3%*%beneficios
trab3=matrix(c(3000,3000,3000,3000,3000,3000,3000),ncol=1)
benef.finmes3=finmes3-trab3
cat("\n Un mes aún no se considera estado estacionario, por lo que si influye por qué máquina se comience, por ello diremos que la empresa generará a fin de mes contando dos sueldos de 1000€ cada uno, un beneficio entre:",max(benef.finmes3),"y",min(benef.finmes3))##
## Un mes aún no se considera estado estacionario, por lo que si influye por qué máquina se comience, por ello diremos que la empresa generará a fin de mes contando dos sueldos de 1000€ cada uno, un beneficio entre: 25907.84 y 23442.56Gráficos representatativos para las conclusiones
x=c(0:6)
plot(estados,benef.finmes1,type="l",col="blue",xlab="Máquinas estropeadas iniciales", ylab="Beneficio Fin de mes")
lines(estados,benef.finmes2,type="l",col="green")
lines(estados,benef.finmes3,type="l",col="red")
legend("bottomleft",col=c("blue","green","red"),legend=c("1 reparador","2 reparadores","3 reparadores"), lwd=3,bty="n")
data=list(benef.finmes1,benef.finmes2,benef.finmes3)
boxplot(data,col=c("blue","green","red"),xlab="Número de reparadores",ylab="Beneficio Final de Mes")
legend("bottomright",col=c("blue","green","red"),legend=c("1 reparador","2 reparadores","3 reparadores"), lwd=3,bty="n")
Conclusiones
Viendo los resultados presentados de forma gráfica un poco más arriba vemos unas conclusiones no tan obvias. Debemos destacar que un mes no representa un estado estacionario, es decir, es muy relevante como podemos ver el número de máquinas estropeadas con el que se comienza.
1 reparador es muy útil cuando el número de máquinas que comienzan estropeadas es igual a 0, pues cuesta que se acumule el trabajo, sin embargo, como vemos en el gráfico, cuando hay más de 1 rotura el beneficio cae drásticamente, pues tiene que ir reparando de uno en uno.
2 reparadores es bajo nuestro punto de vista lo más óptimo para la empresa, pues la curva de beneficio y rotura no sufre de grandes bajadas, y el beneficio se mantiene medianamente estable con cualquier número de reparaciones. En el boxplot podemos ver como este es el que tiene la media de beneficio más alta, por lo que proporciona mucha confianza.
3 reparadores finalmente podemos ver que es un valor inadmisible, pues siempre perderíamos dinero, la tasa de rotura y reparación de las máquinas hace que 3 reparadores sea una pérdida de dinero en cualquier punto comparado con dos reparadores. Aunque a pesar de no ser nada bueno, sí resulta mejor tener 3 reparadores cúando hay 3 máquinas rotas que uno solo. Vemos que 3 reparadores es un claro ejemplo de sobrecostes.
Si conociéramos el porcentaje de veces en las que la producción se inicia en cada estado podríamos ofrecer conclusiones mucho más detalladas, con esto nos referimos a que si supiéramos que por ejemplo el 90% de las veces se inicia con 0 reparaciones deberíamos contratar solo un reparador por ser esto lo más idóneo.
Se han representado también los boxplots para mostrar de forma clara como un trabajador es el valor que más beneficio permite, pero también el que menos, por ser muy eficaz con pocas máquinas en reparación y muy poco con muchas.