Taller 2: RStudio y análisis de falla

Christian Vargas

1. RStudio Scripts Cp y Fallas

Considerando el documento RStudio_Scripts_Cp y Fallas (que se encuentra en la plataforma Classroom del curso en la sección Lenguaje R y RStudio) que contiene Scripts y ejemplos de manejo de datos, replicar los siguientes Scripts:

a. Análisis en R

Análisis estadístico básico en R, interpretando el diagrama boxplot (¿Qué representa cada uno de sus componentes?).

E1

summary(E1)

       X       
 Min.   :2045  
 1st Qu.:2684  
 Median :3456  
 Mean   :3615  
 3rd Qu.:4454  
 Max.   :5890  

boxplot(E1, main="Boxplot de E1", ylab="Valor",
        col="lightblue")

Caja (IQR): del primer al tercer cuartil.

Línea central: mediana.

Bigotes: llegan hasta 1.5·IQR; puntos fuera son atípicos.

b. Capacidad de proceso

library(qcc)

# Datos de diámetros (reemplaza con tus datos de clase)
CP <- c(31.6, 33.4, 35.0, 32.2, 33.8,
       38.4, 31.6, 36.8, 35.0, 34.0,
       29.8, 34.0, 33.0, 34.8, 35.6,
       30.8, 33.0, 31.6, 33.8, 33.8)

# Índice para cada dato
Muestra <- 1:length(CP)

# Histograma
hist(CP, main="Histograma de Diámetros", xlab="Diámetro (mm)", col="lightgreen")

D_groups <- qcc.groups(CP,Muestra)

# Grupos para qcc (aunque no hay grupos reales, usamos los datos tal cual)
q <- qcc(D_groups, type = "xbar", nsigmas = 3, plot = FALSE)

# Capacidad de proceso
process.capability(q, spec.limits = c(30, 38))

Process Capability Analysis

Call:
process.capability(object = q, spec.limits = c(30, 38))

Number of obs = 20           Target = 34
       Center = 33.6            LSL = 30
       StdDev = 2.303           USL = 38

Capability indices:

       Value    2.5%   97.5%
Cp    0.5791  0.3965  0.7614
Cp_l  0.5212  0.3358  0.7066
Cp_u  0.6370  0.4274  0.8465
Cp_k  0.5212  0.3003  0.7421
Cpm   0.5705  0.3926  0.7482

Exp<LSL 5.9%     Obs<LSL 5%
Exp>USL 2.8%     Obs>USL 5%

El proceso tiene una baja capacidad para cumplir consistentemente con las especificaciones establecidas. Esto sugiere que una proporción importante de piezas puede estar fuera de tolerancia, especialmente por debajo del límite inferior.

El índice de capacidad Cp fue de 0.5791, lo que indica que el proceso no es suficientemente capaz, ya que este valor es menor a 1. Además, los índices Cp_l y Cp_u, de 0.5212 y 0.6370 respectivamente, evidencian que el proceso está descentrado respecto a las especificaciones, con una mayor proporción de piezas que tienden a caer fuera del límite inferior (5% de observaciones reales bajo el LSL y 2.8% por encima del USL). Estos resultados indican que el proceso necesita ajustes para mejorar su estabilidad y precisión, que se puede lograr reduciendo la variabilidad.

c. Análisis de fallas en máquinas

Análisis de fallas en máquinas, analizando la condición del equipo considerado.

Horas <- c(2045, 3456, 4567, 3421, 2678, 3789, 4321, 5890, 2367, 4398,
           2690, 4532, 3190, 2367, 4509)

# 3. Análisis estadístico
summary(Horas)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2045    2684    3456    3615    4454    5890 

hist(Horas, main="Histograma de fallas (horas)", xlab="Horas", col="orange", border="white")

# 4. Ajuste a distribución de Weibull
ajuste_weibull <- fitdist(Horas, "weibull")

Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs
Warning :Se han producido NaNs

plot(ajuste_weibull)

# 5. Resumen del ajuste
summary(ajuste_weibull)

Fitting of the distribution ' weibull ' by maximum likelihood 
Parameters : 
Loglikelihood:  -125.5636   AIC:  255.1272   BIC:  256.5433 
Correlation matrix:
          shape     scale
shape 1.0000000 0.3266744
scale 0.3266744 1.0000000

El valor de β indica que el equipo se encuentra en la fase de desgaste, lo que significa que la probabilidad de falla aumenta a medida que transcurre el tiempo de operación. Con esta información, en un entorno real se podría plantear la implementación de planes de mantenimiento preventivo enfocadas al sostenimiento del equipo.

Además, el valor de \(\alpha\)(escala) proporciona una referencia sobre el tiempo promedio en el que se espera la ocurrencia de fallas, lo que permite programar mantenimientos alrededor de las 3000 a 3500 horas, considerando un margen de seguridad para minimizar riesgos de paros no programados

d. Aplicación

Aplicar el script de R con los datos tomados en clase (20 datos) sobre la fabricación del componente mecánico cuyo diámetro fue especificado en \(16,90 \pm 0,08\) mm, Tolerancia \([16,82 – 16,98]\). Analizar la información de los resultados generados y plantear una estrategia de mejoramiento de proceso para lograr la especificación.

# 1. Datos de diámetros
Diametro <- c(16.88, 16.81, 16.92, 16.92, 16.81, 16.93, 16.96, 16.96,
              16.94, 16.95, 16.95, 16.96, 16.92, 16.97, 16.94, 16.94,
              16.93, 16.93, 16.82, 16.93)

# 2. Análisis estadístico
summary(Diametro)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  16.81   16.92   16.93   16.92   16.95   16.97 

boxplot(Diametro, main="Boxplot de Diámetros", col="lightblue", ylab="Diámetro (mm)")

# 3. Histograma
hist(Diametro, main="Histograma de Diámetros", xlab="Diámetro (mm)", col="lightgreen", border="white")

# 4. Capacidad de proceso (usando qcc)
if (!require(qcc)) install.packages("qcc"); library(qcc)

# Control de calidad para valores individuales
q <- qcc(Diametro, type = "xbar.one", nsigmas = 3, labels = 1:20,
         xlab = "Muestra", ylab = "Diámetro (mm)", title = "Carta X-bar de valores individuales")

# Capacidad de proceso (tolerancia 16.82 a 16.98)
process.capability(q, spec.limits = c(16.82, 16.98))

Process Capability Analysis

Call:
process.capability(object = q, spec.limits = c(16.82, 16.98))

Number of obs = 20           Target = 16.9
       Center = 16.92           LSL = 16.82
       StdDev = 0.03873         USL = 16.98

Capability indices:

       Value    2.5%   97.5%
Cp    0.6886  0.4714  0.9054
Cp_l  0.8478  0.5905  1.1051
Cp_u  0.5293  0.3423  0.7164
Cp_k  0.5293  0.3065  0.7522
Cpm   0.6213  0.4137  0.8288

Exp<LSL 0.55%    Obs<LSL 10%
Exp>USL 5.6%     Obs>USL 0%

El resumen estadístico muestra una media cercana a 16,92 mm, lo cual es positivo, pero también resalta la presencia de valores mínimos en 16,81 mm, justo en el límite inferior, y máximos en 16,97 mm, muy cercanos al límite superior. De manera general, se aprecia una buena precisión en el proceso.

Después de realizar los análisis anteriores, se concluye que la capacidad de proceso en algunos casos es baja y el proceso de fabricación puede presentar desviaciones respecto a las especificaciones. Para mejorar, es necesario implementar acciones como la reducción de la variabilidad del proceso, optimización de los parámetros de máquina, ajuste de herramientas, establecimiento de un programa de mantenimiento preventivo, y capacitación del personal en temas de control de calidad

2. Análisis de fallas

A continuación se muestran los datos de tiempo entre fallas TBF, en horas, para un equipo de fabricación.

tbf_horas <- TBF$`TBF (horas)`
cat("TBF (h): ",tbf_horas, "\n")

TBF (h):  14 11 11 17 14 19 16 15 13 15 15 15 13 14 14 15 17 14 14 15 

A partir de estos datos encontrar la función de distribución de Weibull que modela la probabilidad de falla \(F(t)\) y la confiabilidad \(R(t)\).

tbf_fit_weibull <- invisible(rriskFitdist.cont(tbf_horas, "weibull"))
tbf_beta <- invisible(tbf_fit_weibull$estimate["shape"])
tbf_alfa <- invisible(tbf_fit_weibull$estimate["scale"])

tbf_fit_weibull <- invisible(rriskFitdist.cont(tbf_horas, "weibull"))

G2;H2;Warningh in rriskFitdist.cont(tbf_horas, "weibull") :
  Kolmogorov-Smirnov and Anderson-Darling statistics may not be correct with tiesg

tbf_beta <- invisible(tbf_fit_weibull$estimate["shape"])
tbf_alfa <- invisible(tbf_fit_weibull$estimate["scale"])

cat("Escala (alfa):", scale_est, "\nForma (beta):", shape_est, "\n")

Escala (alfa): 15.01125 
Forma (beta): 6.589119 

tbf_dweibull <- invisible(dweibull(tbf_horas, tbf_beta, tbf_alfa))
plot(tbf_horas, tbf_dweibull, col="blue")

Función de densidad de probabilidad

tbf_dens_weibull <- function(x, tbf_beta, tbf_alfa) {
  (tbf_beta / tbf_alfa) * (x / tbf_alfa)^(tbf_beta - 1) * exp(-(x / tbf_alfa)^tbf_beta)
}

A partir de la definición de la función, se procede a graficar para apreciar su comportamiento.

curve(tbf_dens_weibull(x, tbf_beta, tbf_alfa), xlim=c(min(tbf_horas)-0.5,max(tbf_horas)+0.5),xlab="Tiempo entre fallas TBF", ylab="Densidad de Probabilidad", col="red")

¿Cuál será el valor de probabilidad de falla y confiabilidad al cabo de 15 horas de operación? ¿Es adecuado el valor de confiabilidad para la operación del equipo?

---

Probabilidad de falla

tbf_pweibull <- invisible(pweibull(tbf_horas, tbf_beta, tbf_alfa))

plot(tbf_horas, tbf_pweibull, xlab="Tiempo entre fallas TBF",
     ylab="Probabilidad de falla", col="blue")

Se tiene la distribución acumulada de Weibull que dicta la probabilidad de falla:

tbf_acum_weibull <- function(x,tbf_beta,tbf_alfa){
  1 - exp(-(x/tbf_alfa)^tbf_beta)
}

tbf_acum_weibull

function(x,tbf_beta,tbf_alfa){
  1 - exp(-(x/tbf_alfa)^tbf_beta)
}

curve(tbf_acum_weibull(x,tbf_beta,tbf_alfa), xlim=c(min(tbf_horas)-0.5,max(tbf_horas)+0.5))

Dada la curva, si se desea obtener el valor de la función de probabilidad de falla para 15 horas:

prob_falla <- pweibull(15, shape = shape_est, scale = scale_est)

cat("Probabilidad de falla: ",prob_falla*100, "%")

Probabilidad de falla:  63.03032 %

---

Confiabilidad

Dada la probabilidad de falla para 15 horas de operación, la confiabilidad se puede definir como:

confiabilidad <- 1 - prob_falla
cat("Confiabilidad: ",confiabilidad*100, "%")

Confiabilidad:  36.96968 %

Un valor significativamente bajo de confiabilidad como el obtenido indica que el equipo no tiene una capacidad adecuada para desempeñar sus funciones de manera satisfactoria, en la medida que es más propenso a fallar que lo contrario.

Función de tasa de fallas

Expresar matemáticamente la función de tasa de fallas \(\lambda(t)\) y a partir de ésta determinar el valor de la tasa de fallas para 15 horas de operación.`

Recordando que \(\beta =6.59\) y \(\alpha=15.01\), la función lambda de tasa de fallas se presenta entonces como:

\(\lambda(t)=0.53*\left(\frac{t}{15.01}\right)^{5.59}\)

tbf_lambda <- function(x,tbf_alfa,tbf_beta){
  (tbf_beta/tbf_alfa)*(x/tbf_alfa)^(tbf_beta-1)
}

cat("Tasa de fallas para 15 horas", tbf_lambda(15,tbf_alfa,tbf_beta)," fallas/hora")

Tasa de fallas para 15 horas 0.4469101  fallas/hora

Curva de bañera

¿En qué zona de la curva de la bañera se encuentra el equipo? El factor de forma es mayor a 1, por lo que se puede entender que el equipo se encuentra en etapa de envejecimiento. En ese sentido, el equipo se enfrenta a falla por desgaste.

Conclusión

¿Bajo esta condición se consideraría un reemplazo en el corto plazo del equipo?

Entendiendo que

3. Tiempo de reparación (TT)

A continuación, se muestran los datos de los tiempos para reparación o mantenimiento TTR (TTM por algunos autores) para el mismo equipo de la tabla anterior.

ttr_horas <- TTR$`TTR (horas)`
cat("TTR (h): ", ttr_horas,"\n")

TTR (h):  1.7 1.8 4.4 3.3 4.3 3 2.1 2 2.1 1.9 2 1.9 2 2.1 2.1 2.1 1.8 2.2 2.2 2 

Disponibilidad

Determinar la disponibilidad \(A\) (Availability) considerando que el sistema es reparable. ¿Es este valor adecuado para la operación del equipo en la línea de producción?

Se define la disponibilidad, en un sistema reparable, como: \(A=\frac{MTBF}{MTBF+MTTM}\). A continuación se definen estas variables para los conjuntos de datos

MTBF <- mean(tbf_horas)
MTTM <- mean(ttr_horas)

cat("MTBF = ",MTBF," h\nMTTM = ",MTTM," h")

MTBF =  14.55  h
MTTM =  2.35  h

En consecuencia, la disponibilidad A es:

A <- MTBF/(MTBF+MTTM)
cat("A = ",A*100,"%")

A =  86.09467 %

Mantenibilidad

Expresar matemáticamente la función de mantenibilidad para el equipo y determinarla para un TTR de 3 horas.

Asumiendo mantenibilidad con TTM exponencial, se aplicará la ecuación correspondiente:

\(M(t)=1-e^{-\mu t}\)

De donde:

mu <- 1/MTTM
cat("\u03BC = ",mu)

μ =  0.4255319

En consecuencia, la función matemática será:

\(M(t)=1-e^{-0.425t}\)

Y la mantenibilidad M para 3 horas de operación será:

mantenibilidad <- function(x){1-exp(-mu*x)}
cat("M(13) = ", mantenibilidad(3))

M(13) =  0.7210146

curve(mantenibilidad(x), xlim=c(min(ttr_horas)-0.5,max(ttr_horas)+0.5))

Conclusiones

¿Qué se puede concluir sobre la mantenibilidad del equipo?

En cuanto a la mantenibilidad del equipo a partir de los datos planteados se puede apreciar un comportamiento favorable que, aunado a la relativamente alta disponibilidad que dispone, permite un proceso óptimo en tiempos de reparación.