Funciones

# definición de una CMTC en markovchain
# Q matriz generadora
# estados=estados del proceso
# proceso <- new("ctmc", states = estados, byrow = TRUE, 
#             generator = Q, name = "proceso")

# Simulación de ns saltos saliendo de un estado inicial al azar 
# sim=rctmc(n=ns,ctmc=proceso, include.T0=TRUE,out.type="list")

# Simulación hasta Tfin saliendo de un estado inicial dado (o dist.inicial) 
# sim=rctmc(n=Inf,ctmc=proceso, T=Tfin, include.T0=FALSE,
#              out.type="df",initDist)

# Probabilidades de transición (0,T]
# matriz.prob.trans(Rmat, T, cal=1)
# probabilityatT(proceso,T,x0)

# estimación montecarlo
# montecarlo(datosim,type="mean",z=NULL,alpha=0.95)

# tiempos de ocupación
# tiempos.ocupacion<- function(R, T, cal=1)

# comprobamos que una CMTC es irreducible
# is.CTMCirreducible(proceso)

# distribución estacionaria
# steadyStates(proceso)

# Función para la resolución numérica de las ecuaciones de balance
# distr.lim.general<-function(R)

# Tiempos de primer paso
# destino y origen -> numerales
# tiempos.primer.paso(R, destino, estados)
# ExpectedTime(proceso,origen,destino)

Máquinas estropeadas

Para el proceso de Mantenimiento de máquinas en el que se define \(X_t\) como el número de máquinas estropeadas, consideramos cuatro máquinas disponibles y dos operarios para repararlas en caso de fallo, con tiempos de vida de las máquinas exponenciales con media de 3 días (\(Exp(\lambda=1/72)\)), y tiempos medios de reparación de 2 horas \(Exp(\mu=1/2)\). Cuando una máquina se estropea es enviada a reparar, y si todos los reparadores están ocupados, queda en cola de espera hasta que uno se desocupe.

Nota: Todas las tasas están expresadas en horas.

1. Matriz de tasas

Construye la matriz de tasas e interpreta sus elementos. Grafica el sistema.

El espacio de estados es \(S=\{0,1,2,3,4\}\) máquinas estropeadas.

  • Cuando hay 0 máquinas estropeadas, la única opción posible es que alguna de las 4 operativas se averíe. Puesto que los tiempos hasta avería de todas ellas son \(Exp(\lambda)\), el tiempo hasta el primer fallo será el mínimo de todos ellos, que se distribuye \(Exp(4\lambda)\). Así pues, la tasa de avería en este caso, o lo que es lo mismo, de transición al estado 1, es de \(r_{01}=4\lambda\).
  • Cuando hay 1 máquina estropeada, podría ponerse en funcionamiento y pasar el sistema al estado 0 con tasa \(r_{10}=\mu\),o podría estropearse alguna de las otras 3 que están operativas, con tasa \(r_{12}=3\lambda\).
  • Cuando hay 2 máquinas estropeadas, podría ponerse en funcionamiento alguna de ellas con tasa \(r_{21}=2\mu\), o estropearse alguna de las otras dos que están operativas, con tasa \(r_{23}=2\mu\).
  • Cuando hay 3 máquinas estropeadas, podría ponerse en funcionamiento alguna de las dos que están siendo reparadas (pues sólo hay 2 reparadores) con tasa \(r_{32}=2\mu\), o estropearse la otra que queda funcionando, con tasa \(r_{34}=\lambda\).
  • Cuando hay 4 máquinas estropeadas, sólo podría ocurrir que volviera a funcionar alguna de las dos que están siendo reparadas por alguno de los 2 reparadores, con tasa \(r_{43}=2\mu\).

El resto de transiciones no son viables, por lo que su valor es cero. Así la matriz de tasas \(R\) resulta:

# Matriz de tasas
estados <-as.character(0:4)
nestados <- length(estados)

R <- matrix(nrow = nestados, ncol = nestados, data = 0,
            dimnames=list(estados,estados))
lambda <- 1/72
mu <- 1/2 

R[1,2] <- 4*lambda 
R[2,1] <- mu 
R[2,3] <- 3*lambda 
R[3,2] <- 2*mu 
R[3,4] <- 2*lambda 
R[4,3] <- 2*mu 
R[4,5] <- lambda
R[5,4] <- 2*mu
round(R,3)
##     0     1     2     3     4
## 0 0.0 0.056 0.000 0.000 0.000
## 1 0.5 0.000 0.042 0.000 0.000
## 2 0.0 1.000 0.000 0.028 0.000
## 3 0.0 0.000 1.000 0.000 0.014
## 4 0.0 0.000 0.000 1.000 0.000

que deja a la luz que solo son viables las transiciones entre estados colindantes. Se trata pues, de un proceso de nacimiento-muerte, en el que los nacimientos implican averías (incremento del número de máquinas estropeadas) y las muertes representan reparaciones (decrementos del número de máquinas estropeadas).

El diagrama del proceso viene dado por:

RR <- matrix(nrow = nestados, ncol = nestados, data = 0)
RR[1,2] <- "4*lambda"
RR[2,1] <- "mu" 
RR[2,3] <- "3*lambda"
RR[3,2] <- "2*mu"
RR[3,4] <- "2*lambda"
RR[4,3] <- "2*mu" 
RR[4,5] <- "lambda"
RR[5,4] <- "2*mu"
pp <- plotmat(t(RR), pos = 5, curve = 0.5, name = estados,
              lwd = 1, box.lwd = 2, cex.txt = 0.8,
              box.type = "circle", box.prop = 0.5, arr.type = "triangle",
              arr.pos = 0.55, self.cex = 0.6,  
              shadow.size = 0.01, prefix = "", endhead = FALSE, main = "")

2. Tiempos de permanencia

Cuando se estropean dos máquinas, ¿por cuánto tiempo permanecen estropeadas? Cuando están todas funcionando, ¿por cuánto tiempo no se estropea ninguna?

Nos preguntan por los tiempos medios de permanencia en el estado 2 (2 máquinas operativas), \(1/r_2\), y en el estado 0 (todas las máquinas operativas), \(1/r_0\). Estos tiempos medios de permanencia son los inversos a las tasas de permanencia, dado que todos los tiempos en el proceso son exponenciales.

Las tasas de permanencia se calculan sumando las tasas de transición por filas. \[r_i=\sum_{j=1}^{N} r_{ij}.\] Así tendremos:

\[\begin{eqnarray*} r_0 &=& r_{01} \\ r_1 &=& r_{10}+r_{12} \\ r_2 &=& r_{21}+r_{23} \\ r_3 &=& r_{32}+r_{34} \\ r_4 &=& r_{43} \end{eqnarray*}\]

permanencia = apply(R,1,sum)
cat("Tasas de permanencia: ",round(permanencia,3))
## Tasas de permanencia:  0.056 0.542 1.028 1.014 1
tpo.permanencia=1/permanencia

cat("\n Tiempos de permanencia: ",round(tpo.permanencia,3))
## 
##  Tiempos de permanencia:  18 1.846 0.973 0.986 1

Para responder las cuestiones, usamos pues estos tiempos esperados de permanencia. Cada vez que el sistema llega a tener todas las máquinas funcionando (\(X=0\)), continúa así por término medio 18 horas. Cuando llega a tener sólo 2 máquinas operativas, tarda, por término medio 1 horas en cambiar de estado.

3.Probabilidades de salto

Si están estropeadas 3 máquinas, ¿cuál es la probabilidad de que a continuación se estropee una más? ¿Y la de que se arregle una de las dos que están reparando?

Contestamos con las probabilidades de salto \(p_{34}\) (están estropeadas 3 y se estropea una más) y \(p_{32}\) (están estropeadas 3 y se repara 1). Puesto que conocemos las tasas de transición y las de permanencia, calculamos las probabilidades de salto utilizando la relación entre ellas, \[p_{ij}=r_{ij}/r_i\]

psalto <- R/permanencia;
psalto
##           0        1          2          3          4
## 0 0.0000000 1.000000 0.00000000 0.00000000 0.00000000
## 1 0.9230769 0.000000 0.07692308 0.00000000 0.00000000
## 2 0.0000000 0.972973 0.00000000 0.02702703 0.00000000
## 3 0.0000000 0.000000 0.98630137 0.00000000 0.01369863
## 4 0.0000000 0.000000 0.00000000 1.00000000 0.00000000

Así respondemos que la probabilidad de que se estropee otra máquina cuando tenemos 3 máquinas estropeadas es de 0.0136986 y la probabilidad de que una de las que están arreglando vuelva a funcionar es de 0.9863014.

4. Probabilidades de transición

Si la fábrica inició la semana con todas las máquinas operativas:

  1. ¿Cuál es la probabilidad de que transcurridas 8 horas tengamos todas las máquinas en funcionamiento?
  2. ¿Y la probabilidad de que tengamos dos máquinas estropeadas?
  3. ¿Y la de que tengamos a los dos reparadores ocupados?

Respondemos con las probabilidades de transición en \(t=8\) horas que transcurren desde el momento inicial, y en concreto por: ninguna estropeada \(p_{00}(t)\), dos estropeadas \(p_{02}(t)\) y al menos dos estropeadas \(p_{02}+p_{03}+p_{04}\) respectivamente.

# definimos el proceso
estados=as.character(0:4)
Q=R-diag(apply(R,1,sum))
maquinas <- new("ctmc", states = estados,byrow = TRUE, 
               generator = Q,name = "maquinas")
# calculamos las probabilidades
t=8
probabilityatT(maquinas,t)
##           0          1           2            3            4
## 0 0.8978073 0.09806895 0.004014511 0.0001078112 1.411770e-06
## 1 0.8826205 0.11196524 0.005238366 0.0001729243 2.952704e-06
## 2 0.8671344 0.12572079 0.006827022 0.0003095664 8.220433e-06
## 3 0.8383399 0.14940656 0.011144389 0.0010524671 5.670628e-05
## 4 0.7904107 0.18368179 0.021307363 0.0040828523 5.173022e-04
ptran=round(probabilityatT(maquinas,t,x0=1),4)

#1 todas las máquinas funcionando (0 estropeadas)
cat("\n Pr[X(8)=0|X(0)=0]",ptran[1])
## 
##  Pr[X(8)=0|X(0)=0] 0.8978
#2 dos máquinas estropeadas
cat("\n Pr[X(8)=2|X(0)=0]",ptran[3])
## 
##  Pr[X(8)=2|X(0)=0] 0.004
#3 todos los operadores ocupados = al menos 2 máquinas estropeadas 
cat("\n Pr[X(8)>=2|X(0)=0]",sum(ptran[3:5]))
## 
##  Pr[X(8)>=2|X(0)=0] 0.0041

5. Probabilidades límite o estacionarias

Si la fábrica inició la semana con todas las máquinas operativas, ¿cuál es la probabilidad de que al inicio de la semana siguiente tengamos todas las máquinas en funcionamiento? ¿Y al cabo de las dos semanas?

Podemos resolver de nuevo utilizando la matriz de probabilidades de transición, y en particular \(p_{00}(t)\) con \(t=24*7\) y \(t=24*14\) respectivamente, o podemos utilizar la aproximación a largo plazo, esto es, directamente la distribución estacionaria. Comprobamos que los resultados son similares, pues al cabo de una semana la distribución de la cadena ya ha llegado a la estacionariedad; de hecho, todas las filas son similares, pues el proceso ha “olvidado” dónde se inició.

NOTA. Puesto que la convergencia es diferente en cada cadena, es necesario que en cada caso verifiquemos, con la matriz de probabilidades de transición, si realmente se ha producido la convergencia en el tiempo solicitado, esto es, si todas las filas son iguales.

# Probabilidades de transición a una semana vista
t=24*7
probabilityatT(maquinas, t)
##           0          1           2           3            4
## 0 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 1 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 2 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 3 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 4 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
# Distribución estacionaria
steadyStates(maquinas)
##              0          1           2           3            4
## [1,] 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
# Probabilidades de transición a dos semanas vista
t=2*24*7
probabilityatT(maquinas, t)
##           0          1           2           3            4
## 0 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 1 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 2 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 3 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
## 4 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06

6. Tiempos de ocupación

Si la fábrica inició la semana con todas las máquinas operativas y ha estado funcionando durante 8 horas consecutivas: 1. ¿Cuánto tiempo ha tenido todas las máquinas operativas? 2. ¿Cuánto tiempo ha estado ocupado sólo un reparador? 3. ¿Cuánto tiempo han estado ocupados los dos reparadores?

Respondemos con la matriz de los tiempos de ocupación, que obtenemos a través de la función tiempos.ocupacion(), que implementa el método de uniformización para aproximarla mediante una suma finita. Puesto que nos dan el estado inicial, \(X_0=0\) (todas las máquinas funcionando), y el periodo transcurrido, \(t=8\), nos fijaremos en la primera fila de la matriz \(M(t)\), \((m_{00}(t), m_{01}(t),..., m_{04}(t))\).

  1. Todas las máquinas operativas implica ninguna estropeada: \(m_{00}(t)\).

  2. Para saber cuánto tiempo ha estado ocupado sólo un reparador, lo relacionamos con el número de máquinas estropeadas. Esto habrá ocurrido cuando sólo una máquina esté averiada, esto es, \(X=1\), y la solución será \(m_{01}(t)\).

  3. Por último, los dos reparadores habrán estado ocupados cuando hayan tenido 2, 3 o 4 máquinas estropeadas, luego el tiempo total de ocupación será la suma de los tiempos de ocupación de estos estados, esto es, \(m_{02}(t)+m_{03}(t)+m_{04}(t)\).

T=8
tocupa=tiempos.ocupacion(R,T,1)
# comprobamos que cubren las 8 horas
apply(tocupa,1,sum)
##      0      1      2      3      4 
## 8.0001 7.9999 8.0001 8.0000 8.0000
cat("\n Tiempo con todas las máquinas funcionando: ",tocupa[1,1],"horas.")
## 
##  Tiempo con todas las máquinas funcionando:  7.3642 horas.
cat("\n Tiempo con sólo un reparador ocupado: ",tocupa[1,2],"horas.")
## 
##  Tiempo con sólo un reparador ocupado:  0.6139 horas.
cat("\n Tiempo con los dos reparadores ocupados: ",sum(tocupa[1,3:5]),"horas.")
## 
##  Tiempo con los dos reparadores ocupados:  0.022 horas.

7. Proporción del tiempo de ocupación

Si la fábrica inició la semana con todas las máquinas operativas y ha estado funcionando durante 8 horas consecutivas,

  1. ¿qué fracción del tiempo ha tenido todas las máquinas operativas?
  2. ¿Cuál es el aprovechamiento total de los dos reparadores durante una jornada laboral? ¿Cuántas horas han trabajado en relación a las jornadas laborales que tienen ambos?
  3. ¿Cuánto sueldo podríamos haber ahorrado?

Utilizamos los resultados que obtuvimos en el apartado anterior, y respondemos en términos de fracciones de tiempo (dividiendo por el tiempo de funcionamiento del sistema \(T=8\)).

Para responder la pregunta 2, el tiempo de trabajo total dedicado a reparaciones por los dos reparadores será igual al tiempo en que la fábrica sólo ha tenido una máquina averiada, periodo en el que ha trabajado sólo un reparador, más el tiempo en que la fábrica ha tenido dos o más máquinas averiadas, periodo en el cual han trabajado los dos reparadores. El tiempo total de trabajo de los dos trabajadores es \(T=8 \cdot 2\) horas, juntando sus dos jornadas laborales. El tiempo total de trabajo efectivo de los dos trabajadores será \[t=m_{01}(8)+2 \cdot [m_{02}(8)+m_{03}(8)+m_{04}(8)]\] La proporción de horas de trabajo que han trabajado será pues \(t/(2\cdot 8)\).

#1 todas las máquinas operativas
T=8
cat("\n Tiempo con todas las máquinas funcionando: ",tocupa[1,1]/T*100,"%")
## 
##  Tiempo con todas las máquinas funcionando:  92.0525 %
#2  tiempo de trabajo total de los dos reparadores
cat("\n Tiempo de trabajo realizado en 8 horas: ",sum(tocupa[1,2:5])+sum(tocupa[1,3:5])*2,"horas")
## 
##  Tiempo de trabajo realizado en 8 horas:  0.6799 horas
trabajo=(sum(tocupa[1,2:5])+sum(tocupa[1,3:5])*2)/(2*T)*100
cat("\n Aprovechamiento de las dos jornadas laborales: ",trabajo,"%")
## 
##  Aprovechamiento de las dos jornadas laborales:  4.249375 %
cat("\n % Ahorrable de las dos jornadas laborales: ",100-trabajo,"%")
## 
##  % Ahorrable de las dos jornadas laborales:  95.75062 %

8. Valor resperado

¿Cuántas máquinas esperamos que estén operativas cuando hayan transcurrido 8 horas desde el inicio de la jornada, si se empezó con todas las máquinas operativas?

Puesto que \(X_8\) representa el número de máquinas estropeadas después de 8 horas, \(4-X_8\) será el número de máquinas funcionando después de 8 horas. Nos preguntan pues, por \(4-E[X_8|X_0=0]\), valor esperado que podemos calcular con las probabilidades de transición \(p_{0j}(8)\), esto es, \[E[X_8|X_0=0]=\sum_{j=0}^4 j \cdot p_{0j}(8).\]

T=8
ptran <- probabilityatT(maquinas,T,1)
estados<-0:4
estro8=ptran%*%estados; estro8
##          [,1]
## [1,] 0.106427
cat("Núm.Esperado de máquinas funcionando tras 8 horas:",4-estro8)
## Núm.Esperado de máquinas funcionando tras 8 horas: 3.893573

9. Estado estacionario

** Si la fábrica inició la semana con todas las máquinas operativas y ha estado funcionando durante una semana entera:

  1. ¿cuánto tiempo ha tenido todas las máquinas operativas?
  2. ¿Qué fracción del tiempo ha tenido todas las máquinas operativas?
  3. En general, ¿qué porcentaje del tiempo de trabajo en la fábrica están todas las máquinas operativas?

Podemos responder a esta pregunta reproduciendo el resultado en apartados previos, con los tiempos de ocupación, y en particular con \(m_{00}(t)\) para \(t=24 \cdot 7\). Sin embargo, puesto que una semana se considera mucho tiempo, el comportamiento que observemos se parecerá mucho al comportamiento a largo plazo del proceso, es decir, los tiempos que nos piden provendrán de la “distribución de ocupación del proceso”, que coincide con la distribución límite y con la distribución estacionaria \(p=(p_0,p_1,...,p_4)\), y en particular de \(p_0\). Comprobamos no obstante, que los resultados son similares.

# tiempos de ocupación
t=24*7
tocupa=tiempos.ocupacion(R,t)[1,];tocupa
##        0        1        2        3        4 
## 150.7531  16.5427   0.6850   0.0189   0.0003
cat("\n Tiempo con todas las máquinas funcionando en una semana: ",round(tocupa[1],2),"=",round(tocupa[1]/t*100,2),"%")
## 
##  Tiempo con todas las máquinas funcionando en una semana:  150.75 = 89.73 %
# distribución límite
p.estacionaria=round(steadyStates(maquinas),4)
cat("\n Tiempo con todas las máquinas funcionando en cualquier periodo a largo plazo: ",round(p.estacionaria[1]*100,2),"%")
## 
##  Tiempo con todas las máquinas funcionando en cualquier periodo a largo plazo:  89.62 %

10. Análisis de costes con inicio dado

Si sabemos que cada minuto que una máquina está averiada genera unos costes de 1€ y cuando está operativa genera unos beneficios de 5€/hora, ¿cuánto gana la empresa en una jornada de 8 horas, si esta empieza con todas las máquinas funcionando? (No consideramos el sueldo de los reparadores (y por lo tanto el coste de cada reparación).

Los beneficios esperados acumulados durante un periodo de \(T=8\) horas, sin considerar el sueldo de los reparadorees, los cuantificamos con la fórmula \[g(T)=m(T) \cdot c,\] donde \(m(T)\) es la matriz de ocupación y \(c\) una matriz columna con los beneficios incurridos en cada estado por unidad de tiempo, por avería y funcionamiento. Obtenemos un vector con los beneficios obtenidos en función de cuál es el estado inicial del sistema, \((g(0,T), g(1,T),g(2,T),g(3,T),g(4,T))\).

coste.hora=1*60
beneficio.hora=5
beneficios=(4:0)*beneficio.hora-(0:4)*coste.hora; beneficios
## [1]   20  -45 -110 -175 -240
t=8
g=tiempos.ocupacion(R,t,1) %*% matrix(beneficios,ncol=1);g
##        [,1]
## 0  117.2060
## 1   -7.2925
## 2 -133.3885
## 3 -320.3890
## 4 -564.7500
cat("Beneficio esperado de una jornada de 8 horas empezando con todas las máquinas operativas (sin reparadores):",round(g[1,],2),"€")
## Beneficio esperado de una jornada de 8 horas empezando con todas las máquinas operativas (sin reparadores): 117.21 €

11. Análisis de costes con incertidumbre en el inicio

Tras un estudio realizado se ha comprobado que el 89.7% de las semanas se inician con todas las máquinas operativas, el 9.88% con una máquina estropeada, el 0.41% con 2 máquinas estropeadas, el 0.01% con 3 máquinas estropeadas y el 0% con las cuatro máquinas estropeadas. ¿Cuánto gana la empresa en una jornada cualquiera de 8 horas?

Nos preguntan por el beneficio esperado acumulado en 8 horas, considerando la incertidumbre existente sobre cómo se iniciará una jornada cualquiera de 8 horas, y representada por la distribución inicial \(p_0(0),p_1(0),..., p_4(0)\), que se da en el enunciado. Así la respuesta la obtenemos con: \[E[g(t)]=\sum_{j=0}^4 g(j,t) p_j(0)\]

pini=c(89.7,9.88,0.41,0.01,0)/100
coste.medio=as.vector(g) %*% pini
cat("Beneficio esperado total de una jornada cualquiera de 8 horas (sin reparadores):",round(coste.medio,2),"€")
## Beneficio esperado total de una jornada cualquiera de 8 horas (sin reparadores): 103.83 €

12. Análisis de costes a largo plazo

Con los mismos costes especificados previamente para los tiempos de avería y de funcionamiento, y considerando que el sueldo de los reparadores es de 1200€/mes por reparador:

  1. ¿Cuál será la ganancia esperada por producción (por hora) para la empresa, sin considerar el sueldo de los reparadores?
  2. ¿Cuál será la ganancia esperada (por hora) para la empresa, considerando el sueldo de los dos reparadores?
  3. ¿Qué pasaría si los costes por unidad de tiempo que está averiada una máquina se duplicaran respecto de la situación actual?
  4. ¿Qué sueldos recomendarías a la empresa para garantizar unos beneficios de al menos el 70% de las ganancias que se generan por la producción?
  5. Si se decidiera subcontratar las reparaciones asumiendo un contrato de x horas de trabajo (y un máximo de dos reparaciones simultáneas), ¿cuántas horas (máximas) habría de contratar? ¿A cuánto debería pagar la hora para garantizar unos beneficios netos del 95% de los beneficios que genera la producción?

Para responder la primera pregunta multiplicamos los beneficios por la distribución estacionaria, pues es la que regula el funcionamiento del sistema a largo plazo.

coste.hora=1*60
beneficio.hora=5
beneficios=(4:0)*beneficio.hora-(0:4)*coste.hora
g=steadyStates(maquinas)%*%beneficios
cat("Beneficio hora (sin reparadores):",round(g,2),"€")
## Beneficio hora (sin reparadores): 12.97 €
cat("\n Beneficio mes (sin reparadores):",round(g*30*24,2),"€")
## 
##  Beneficio mes (sin reparadores): 9335.15 €

Pagar dos sueldos mensuales de 2000€, implica detraer esos sueldos a las ganancias conseguidas, esto es,

sueldo=1200
cat("\n Beneficio mes (con reparadores):",round(g*30*24-sueldo*2,2),"€")
## 
##  Beneficio mes (con reparadores): 6935.15 €

Si los costes por tiempo averiado se duplicaran, las ganancias de la empresa serían

coste.hora=2*60
beneficio.hora=5
beneficios=(4:0)*beneficio.hora-(0:4)*coste.hora
g=steadyStates(maquinas)%*%beneficios
cat("Beneficio hora (sin reparadores):",round(g,2),"€")
## Beneficio hora (sin reparadores): 6.47 €
cat("\n Beneficio mes (sin reparadores):",round(g*30*24,2),"€")
## 
##  Beneficio mes (sin reparadores): 4659.9 €
cat("\n Beneficio mes (con reparadores):",round(g*30*24-sueldo*2,2),"€")
## 
##  Beneficio mes (con reparadores): 2259.9 €

Así, si queremos garantizar a la empresa unas ganancias al menos del 70% de las ganancias que genera la producción, solo estaría dispuesto a pagar en sueldos al reparadores el 30% restante, por lo que repartiendo entre dos, tendríamos:

ganancias.produccion.mes = g*30*24
cat("\n Beneficio mes (sin reparadores):",round(ganancias.produccion.mes,2),"€")
## 
##  Beneficio mes (sin reparadores): 4659.9 €
sueldo.reparador=0.3*ganancias.produccion.mes/2
cat("\n Sueldo recomendado operador/mes:",round(sueldo.reparador,2))
## 
##  Sueldo recomendado operador/mes: 698.98

13. Tiempos de primer paso

1. Cuando hay dos máquinas estropeadas, ¿cuánto tiempo pasa hasta que están todas estropeadas? 2. Si la fábrica inicia con todas las máquinas operativas, ¿cuánto tiempo pasará hasta que hayan de trabajar simultáneamente los dos operadores?

Respondemos con los tiempos de primer paso desde un estado a otro.

Para que trabajen simultáneamente los dos operadores, el número de máquinas estropeadas habrá de ser al menos de 2, esto es, 2, 3 o 4. Calculemos pues, los tiempos de primer paso desde el estado 0 hasta los estados 2, 3 y 4.

#1 
et=ExpectedTime(maquinas,3,5)
cat("Con dos máquinas estropeadas, trancurren",et,"horas hasta que hay 4 máquinas averiadas, esto es,",round(et/24,1),"días.\n")
## Con dos máquinas estropeadas, trancurren 633420 horas hasta que hay 4 máquinas averiadas, esto es, 26392.5 días.
A=3:5 # corresponde a las filas que identifican los estados 2, 3,4
estados=0:4
primer.paso=tiempos.primer.paso(R,A,estados)
cat("Si iniciamos con todas las máquinas operativas, trancurren",round(primer.paso[1,],1),"horas hasta que los dos reparadores trabajan a la vez, esto es,",round(primer.paso[1,]/24,1),"días.")
## Si iniciamos con todas las máquinas operativas, trancurren 258 horas hasta que los dos reparadores trabajan a la vez, esto es, 10.8 días.

14. Simulación

Queremos aproximar por MC (y comparar con los valores teóricos):

Después de una semana:

  1. ¿Cuántas máquinas estarán en funcionamiento al inicio de la semana siguiente?
  2. ¿Cuánto tiempo están funcionando a la vez las cuatro máquinas (en términos porcentuales)?
  3. Si con 2 máquinas estropeadas la fábrica ha de parar la producción, ¿qué proporción del tiempo para la fábrica durante una semana?
  4. ¿Cuánto tiempo se ha trabajado reparando averías? ¿Cuántas horas de trabajo retribuído no se ha trabajado (en términos porcentuales)?
  5. ¿Cuántas horas están trabajando simultáneamente los dos reparadores (en términos porcentuales durante una jornada laboral)?
  6. ¿Qué recomendación harías sobre cómo gestionar las reparaciones y los pagos a los reparadores?
# definición del proceso
mantenimiento=function(N,M,lambda,mu){
  # N=nº máquinas
  # M=nº reparadores
estados_num=0:N
estados=as.character(estados_num)
R <- matrix(nrow = nestados, ncol = nestados, data = 0,
            dimnames=list(estados,estados))
lambdas=lambda*(N-estados_num[1:N])
min_M=function(x){
  min(x,M)
}
mus=mu*sapply(estados_num[2:(N+1)],min_M)
for(i in 1:N){
  R[i,i+1]=lambdas[i]
  R[i+1,i]=mus[i]
}
Q=R-diag(apply(R,1,sum))
maquinas <- new("ctmc", states = estados,byrow = TRUE, 
               generator = Q,name = "maquinas")
return(maquinas)
}
# Especificación
N=4
M=2
lambda <- 1/72
mu <- 1/2 

maquinas=mantenimiento(N,M,lambda,mu)

nsim=5000
Tsim=7*24
# nº máquinas estropeadas al final del periodo
maq_estropeadas=c()
# tpo funcionamiento simultáneo de todas las máquinas
tpo_fun_full=c()
# tpo en que para la fábrica (X>=2)
tpo_paro=c()
# tpo dedicado a reparaciones
tpo_repara=c()
# tpo de trabajo simultáneo de todos los reparadores
tpo_simul=c()
# función para calcular el mínimo de un número y el num.reparadores
min_M=function(x){
  min(x,M)}

for(i in 1:nsim){
  sim=rctmc(n=Inf,ctmc=maquinas,T=Tsim,out.type = "df",
          initDist=c(1,rep(0,N)),include.T0=TRUE)
  #estados
  sim$states=as.numeric(sim$states)
  # instantes de transición
  sim$time=as.numeric(sim$time)
  m=nrow(sim)
  # tiempos de permanencia
  if(m==1){
    sim$stay=Tsim-sim$time
  }else{
     sim$stay=c(sim$time[2:m],Tsim)-sim$time[1:m]
  }
  
 # sumarios de interés
  maq_estropeadas[i]=tail(sim$states,1)
  tpo_fun_full[i]=sum(sim$stay[sim$states==0])
  tpo_paro[i]=sum(sim$stay[sim$states>=2])
  
  sim$maq_repara=sapply(sim$states,min_M)
  tpo_repara[i]=sum(sim$maq_repara*sim$stay)
  
  tpo_simul[i]=sum(sim$stay[sim$states>=M])
}
## 
##  Nº máquinas estropeadas al iniciar la semana siguiente
## Estimación MC: 0.1118 [ 0.1027714 , 0.1208286 ]
## 
##  % tiempo que funcionan todas las máquinas a la vez
## Estimación MC: 89.74943 [ 89.62502 , 89.87385 ]
## 
##  % tiempo que para la fábrica
## Estimación MC: 0.4076085 [ 0.3879018 , 0.4273153 ]
## 
##  Tiempo trabajado en reparaciones
## Estimación MC: 17.90573 [ 17.68235 , 18.12912 ]
## 
##  % Tiempo retribuido sin trabajar
## Estimación MC: 94.67091 [ 94.60443 , 94.7374 ]
## 
##  % Tiempo que trabajan a la vez todos los reparadores
## Estimación MC: 0.4076085 [ 0.3879018 , 0.4273153 ]

Máquinas operativas

Para el proceso de Mantenimiento de máquinas en el que se define \(X_t\) como el número de máquinas en funcionamiento, consideramos cuatro máquinas disponibles y dos operarios para repararlas en caso de fallo, con tiempos de vida de las máquinas exponenciales con media de 3 días (\(Exp(\lambda=1/72)\)), y tiempos medios de reparación de 2 horas \(Exp(\mu=1/2)\).

  1. Estamos interesados en calcular el valor esperado del número de máquinas en funcionamiento después de 8 horas de funcionamiento, suponiendo que el sistema empieza con todas las máquinas funcionando.
  2. La probabilidad a largo plazo de que todas las máquinas estén en funcionamiento,

Así la matriz de tasas \(R\) resulta:

# Matriz de tasas
estados <- as.character(0:4)
nestados <- length(estados)

R_fun <- matrix(nrow = nestados, ncol = nestados, data = 0,
            dimnames=list(estados,estados))
mu <- 1/2 # tasa reparación
lambda <- 1/72  # tasa funcionamiento

R_fun[1,2] <- 2*mu
R_fun[2,1] <- lambda 
R_fun[2,3] <- 2*mu
R_fun[3,2] <- 2*lambda
R_fun[3,4] <- 2*mu
R_fun[4,3] <- 3*lambda 
R_fun[4,5] <- mu
R_fun[5,4] <- 4*lambda
RR <- matrix(nrow = nestados, ncol = nestados, data = 0)
RR[1,2] <- "2*mu"
RR[2,1] <- "lambda" 
RR[2,3] <- "2*mu"
RR[3,2] <- "2*lambda"
RR[3,4] <- "2*mu"
RR[4,3] <- "3*lambda" 
RR[4,5] <- "mu"
RR[5,4] <- "4*lambda"
pp <- plotmat(t(RR),  pos = 5, curve = 0.5, name = estados,
              lwd = 1, box.lwd = 2, cex.txt = 0.8,
              box.type = "circle", box.prop = 0.5, arr.type = "triangle",
              arr.pos = 0.55, self.cex = 0.6,  
              shadow.size = 0.01, prefix = "", endhead = FALSE, main = "")

1. Valor esperado

Para calcular el valor esperado del número de máquinas en funcionamiento tras 9 horas de funcionamiento,calculamos la matriz de probabilidades de transición para tener la distribución de transición en 9 horas arrancando de 4 máquinas funcionando.

# definimos el proceso
estados=as.character(0:4)
Q_fun=R_fun-diag(apply(R_fun,1,sum))
maquinas_fun <- new("ctmc", states = estados,byrow = TRUE, 
               generator = Q,name = "maquinas_fun")
# calculamos las probabilidades
t=8
estados_num=0:4
prob9_from4=probabilityatT(maquinas,t)[5,]
# valor esperado
esperanza <- round(estados_num%*%prob9_from4, 1)
cat("E[X(9)]=Nº esperado de máquinas funcionando tras 9 horas=",round(esperanza,2))
## E[X(9)]=Nº esperado de máquinas funcionando tras 9 horas= 0.2

2. Distribución estacionaria

estacionaria=steadyStates(maquinas_fun)
estacionaria
##              0          1           2           3            4
## [1,] 0.8961608 0.09957343 0.004148893 0.000115247 1.600653e-06
cat("p_4=Probabilidad todas las máquinas operativas=",round(estacionaria[5],5))
## p_4=Probabilidad todas las máquinas operativas= 0
cat("p_0=Probabilidad todas las máquinas estropeadas=",round(estacionaria[1],5))
## p_0=Probabilidad todas las máquinas estropeadas= 0.89616