Documentación acerca de estimación de

parámetros para TBF en RMES y posibles soluciones para el problema de

elección de mejor distribución explicativa del indicador Software RMES

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Viña del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 empresa@mes.cl Enero, 2011

ÍNDICE  RMES

Índice

1. Análisis Weibull 3 2. Análisis Exponencial 9 3. Soluciones Elección Mejor Distribución 10 4. Apéndice 11

4.1. Aproximación Ecuación 6 . . . . . . . . . . . . . . . . . . . . . . 11 4.2. Regresión Lineal Simple . . . . . . . . . . . . . . . . . . . . . . . 12

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 2

1 ANÁLISIS WEIBULL  RMES

1. Análisis Weibull

Dada cualquier función g(x) con g(0) = 0 y crezca monótonamente hacia el infinito mientras x tienda al infinito, es posible definir una función de probabilidad acumulada como:

F (t) = 1 − e^{−g(t)}

Obviamente, la probabilidad en el tiempo t=0 es cero y crece monótonamente hasta 1 mientras t tienda al infinito. La correspondiente función de densidad f(t) es la derivada de F (t):

f(t) = g′(x)e^{−g(t)}

La tasa de falla para una función de densidad dada se define como

R(t) = \frac{f(t)}{1 − F (t)}

entonces, la tasa de falla para la función de densidad anterior es:

R(t) = \frac{g′(x)e^{−g(t)}}{1 − (1 − e^{−g(t)})} = g′(x)

Esto nos permite definir una distribución con una tasa de ocurrencia R(t) especificada. Una familia de funciones de tasa de 2 parámetros muy útil es:

R(t) = \frac{α}{β}(\frac{t}{α})^{β−1}  (1)

donde α y β son constantes. Esta familia de distribuciones, llamadas así por el ingeniero Suizo Waloddi Weibull (1887-1979) quién popularizó su uso para análisis de confiabilidad, especialmente para modos de fallas metalúrgicos. El primer paper de Weibull acerca del tema fue publicada en 1939, pero el método no atrajo mucho hasta el año 1950. La ”distribución Weibull” también fue estudiada en los años 20 por el estadísta Emil Gumbel (1891-1966), que es recordado hoy por su enfrentamiento con los nazis en 1931 cuando se organizó una campaña para obligarlo a salir de su cátedra en la Universidad de Heidelberg por su abierto opinión pacifista y puntos de vista anti-nazi.

La constante α es llamada parametro de escala, porque escala a la variable t, y la constante β es llamada parámetro de forma, porque determina la forma de la función de tasa. (Ocasionalmente, la variable t en la definición anterior es reemplazada por t-g, donde g es un tercer parámetro usado para definir un punto cero adecuado). Respecto del parámetro de forma β:

Si β > 1 la tasa crece con t

Si β < 1 la tasa decrece con t

Si β = 1, la tasa es constante.

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 3

Figura 1: Tasa de falla

En el último caso, la distribución Weibull se iguala a una distribución Exponencial. La forma de la función tasa de falla se ilustra en la figura 1.

Dado que R(t) es igual a g′(t), integramos esta función para obtener:

g(t) = (\frac{t}{α})^{β}

Claramente, para cualquier valor positivo de α y β la función g(t) crece monótonamente desde 0 hasta el infinito mientras t se lo haga también, lo que nos lleva a definir una probabilidad de densidad acumulada:

F (t) = 1 − e^{−( \frac{t}{α})^{β}}  (2)

y su correspondiente función de densidad:

f(t) = \frac{β}{α}(\frac{t}{α})^{β−1}e^{−( \frac{t}{α})^{β}}  (3)

Supongamos que tenemos un conjunto de n procesos, donde n es un número grande, cada uno comienza operando continuamente al mismo tiempo t = 0 (entiéndase proceso como intervalo de tiempo de buen funcionamiento; con esto se pretende simplificar el problema del funcionamiento continuo de un equipo, asumiendo que cada vez que éste falle y reestablezca su funcionamiento la variable de tiempo t vuelve a cero, de esta forma t no representa tiempo de continuidad asociada a fecha, si no que a la duración de los procesos). Cada proceso tiene una distribucion acumulada de falla de tipo Weibull dada por la ecuación 2 para parámetros fijos α y β, entonces

El numero esperado N(t) de procesos para el tiempo t (es decir, del conjunto n de procesos, la cantidad de ellos que aún seguirían en funcionamiento hasta el tiempo t) es:

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 4

N(t) = F (t)n  (4)
N(t) = (1 − e^{−( \frac{t}{α})^{β}})n  (5)

Reordenando la ecuación 5:

1 − \frac{N(t)}{n} = e^{−( \frac{t}{α})^{β}}

Aplicando ln(.) y multiplicando por −1

ln (\frac{1}{1 − \frac{N(t)}{n}} ) = (\frac{t}{α})^{β}

Aplicando nuevamente ln(.):

ln ( ln (\frac{1}{1 − \frac{N(t)}{n}} )) = β ln(t) − β ln(α)  (6)

Una aproximación a la ecuacion 6, desarrollada en la sección 4.1 del Apéndice, es:

ln (\frac{N(t)}{n}) = β ln(t) + β ln(α)  (7)

Sin embargo, es más fácil usar la ecuación exacta 6 en vez de la ecuación aproximada 7

Dada un conjunto inicial de n = 100 de procesos que comiencen en t = 0 y acumulando horas continuamente posteriormente, supongamos que la primera falla ocurre en un tiempo t = t_{1}. Aproximadamente, podríamos decir que el número esperado de fallas al tiempo de la primera falla es 1 (o sea, aprox N(t_{1}) = 1, entonces, según la ecuación 4, F (t_{1}) = N(t_{1})/n = 1/100. Sin embargo, esto (decir que N(t_{1}) = 1) no es muy óptimo, dado que estadísticamente la primera falla probablemente ocurra algo antes de que el número esperado de fallas alcance 1. Para enteder el porqué de ésto, consideremos un conjunto que consista únicametne de un sólo proceso, en tal caso, el número esperado de fallas a un tiempo t cualquiera sería meramente F (t), el cual sólo se aproxima a 1 mientras t tienda al infinito, y sin embargo, el tiempo medio de falla es en el tiempo t = t_{mediana} donde F (t_{mediana}) = 0,5. En otras palabras, hay una probabilidad de 0.5 de que la falla ocurra antes que el tiempo medio y de un 0.5 de que ocurra después. De ahí que para un conjunto de tamaño 1, el número esperado de fallas al tiempo medio de la primera falla es sólo de 0.5.

En general, dado un conjunto de n procesos, cada uno con una función de densidad f(t), la probabilidad individual de que alguno falle al tiempo t_{m} es F (t_{m}) = N(t_{m})/n. Definiendo este valor por Φ, la probabilidad exacta de que j procesos fallen y n − j no lo hagan al tiempo t_{m} sigue una distribución binomial [4] y es:

P (j; n) = \binom{n}{x}Φ^{j}(1 − Φ)^{n−j}

Y por consiguiente, la probabilidad de que j o más fallen al tiempo t_{m} es:

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 5

P (≥ j; n) = ∑_{i=j}^{n} \binom{n}{i}Φ^{j}(1 − Φ)^{n−i}

Esto representa la probabilidad de que la j-ésimo (de n) falla ocurra en el tiempo t_{m}, y por supuesto, el complemento es la probabilidad de que la j-ésimo aún no ocurra al tiempo t_{m}.

Por lo tanto, dado que la j-ésima falla ocurre al tiempo t_{m}, el valor ”medio” de F (t_{m}) = Φ se obtiene igualando P (≥ j; n) = 0,5 en la ecuación de arriba y resolviendo para Φ. Este valor se conoce como el rango medio y puede ser computado numéricamente. Un enfoque alternativo es utilizar la muy buena aproximación:

F (t_{m}) = \frac{j − 0,3}{n + 0,4}  (8)

Este es el valor (en vez de j/n) que se debería asignar a N(t_{j})/n para la j-ésima falla.

Luego, podemos reemplazar este valor en la ecuación 6, obteniendo

ln ( ln (\frac{1}{1 − \frac{j−0,3}{n+0,4}} )) = β ln(t) − β ln(α)  (9)

Utilizando regresión simple, detallada en la sección del Apéndice 4.2, y usando como variables auxiliares

x_{j} = ln(t_{j})

y_{j} = ln ( ln (\frac{1}{1 − \frac{j−0,3}{n+0,4}} ))

Obtenemos

y_{j} = βx_{j} − β ln(α)

Luego

β = \frac{k ∑ x_{j}y_{j} − ∑ x_{j} ∑ y_{j}}{k ∑ x_{j}^{2} − (∑ x_{j})^{2}}  (10)

α = exp[\frac{∑ y_{j} ∑ x_{j}^{2} − ∑ x_{j} ∑ x_{j}y_{j}}{−β(k ∑ x_{j}^{2} − ( ∑ x_{j})^{2})} ]  (11)

donde k es la cantidad de datos.

La dependencia de la tasa de falla R(t) en el tiempo de cada proceso es también una razón que consideramos en el conjunto coherente de procesos cuyas edades están sincronizadas (t_{0} = 0). Esto simplifica enormemente el análisis. En situaciones más realistas, el conjunto de proceso cambia, y la ”edad” de cada proceso en el conjunto será diferente, de igual manera, la tasa con la cual las horas operacionales se acumulan es una función de tiempo. Mas generalmente, podríamos considerar un conjunto de procesos tales que cada uno posee su ”tiempo propio”:

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 6

t_{j} = m_{j}(t − g_{j})

para todo t mas grande que g_{j}, donde t es el tiempo, g_{j} es la fecha de nacimiento del j-ésimo proceso y m_{j} es factor de uso operacional. Este tiempo propio es entonces la variable de tiempo para la función de densidad de la Weibull para el j-ésimo proceso y la tasa de fallas genera para todo el conjunto en un tiempo dado se compone de todas las tasas individuales. En un conjunto no coherente, cada proceso tiene distintas distribuciones de falla.

Figura 2: Procesos coherentes con censura a la derecha

La figura 2 muestra procesos coherentes. Los procesos 2, 4 y 5 representan procesos detenidos debido a fallas. El resto, continuan su funcionamiento acumulando horas de funcionamiento a sus respectivas tasas. Ellos son llamados puntos ”censurados a la derecha” o ”suspendidos”, dado que imaginamos que su funcionamiento fue suspendido antes de que el proceso fallara. No sabemos cuando realmente podrían haber fallado, por lo que no pueden ser usados directamente como puntos para el ajuste de la distribución, pero de todas maneras podemos utilizar el hecho de que ellos acumulan horas sin fallas. La aproximación usual en análisis de confiabilidad es primero rankear todos los datos de acuerdo a sus horas acumuladas, tal como muestra la figura 3.

Figura 3: Rankeo de datos para ajuste Weibull

A continuación, asignamos un ajuste al ranking para los procesos que efectivamente falla. Si denotamos k_{j} al ranking respecto a todos los datos y r(j) al ranking ajustado asociados a los procesos que fallan (definiendo r(0) = 0), el ranking ajustado de la j-ésima falla viene dado por la siguiente fórmula:

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 7

r(j) = r(j − 1) + \frac{N + 1 − r(j − 1)}{N + 1 − (k_{j} − 1)}

Luego, para el ejemplo de la figura 2 tenemos:

r(1) = r(0) + \frac{10 + 1 + r(0)}{10 + 1 − (6 − 1)} = 1,667

r(2) = r(1) + \frac{10 + 1 + r(1)}{10 + 1 − (7 − 1)} = 3,533

r(3) = r(2) + \frac{10 + 1 + r(2)}{10 + 1 − (10 − 1)} = 7,267

Y estos son los valores que debemos usar para el ajuste Weibull. Luego, en la ecuacion 12 en vez de utilizar j, utilizamos r(j) asociado a los procesos que efectivamente fallaron.

F (t_{m}) = \frac{r(j) − 0,3}{n + 0,4}  (12)

Finalmente, nuestras variables auxiliares para al ajuste se modifican.

x_{j} = ln(t_{j})

y_{j} = ln ( ln (\frac{1}{1 − \frac{r(j)−0,3}{n+0,4}} ))

Dejando la ecuacion lineal 9 como

ln ( ln (\frac{1}{1 − \frac{r(j)−0,3}{n+0,4}} )) = β ln(t) − β ln(α)  (13)

Y las ecuaciones 11 y 10 para los parámetros α y β respectivamente no se modifican.

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 8

2 ANÁLISIS EXPONENCIAL  RMES

2. Análisis Exponencial

El análisis para esta distribución es símil a la distribución Weibull pero forma mucho más simplificada dado que esta distribución es un caso especial de la Weibull cuando β = 1. Con ello, la tasa de falla se define como:

R(t) = \frac{1}{λ}  (14)

la cual es una tasa constante. Luego, la función de probabilidad acumulada se define como:

F (t) = 1 − e^{−λt}  (15)

y su correspondiente función de densidad es:

f(t) = λe^{−λt}  (16)

Notar que las fórmulas análogas con Weibull simplemente se reemplaza β = 1 y α = 1/λ.

Para el ajuste, realizamos los mismos procedimientos que en el análisis Weibull:

N(t) = (1 − e^{−λt})n

1 − \frac{N(t)}{n} = e^{−λt}

ln (1 − \frac{N(t)}{n}) = −λt

Utilizando el mismo concepto de rankeo explicado anteriormente obtenemos:

ln (1 − \frac{r(j) − 0,3}{n + 0,4}) = −λt  (17)

Utilizando las siguientes variables auxiliares:

x = t

y = ln (1 − \frac{r(j) − 0,3}{n + 0,4})

Obtenemos simplemente:

y = −λx

Y aplicando regresión lineal, detallada en la sección del Apéndice 4.2, obtenemos:

λ = −\frac{∑ xy}{∑ x^{2}}  (18)

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 9

3 SOLUCIONES ELECCIÓN MEJOR DISTRIBUCIÓN  RMES

3. Soluciones Elección Mejor Distribución

Lamentablemente, los métodos de estimación explicados en las secciones anteriores no permiten utilizar un Test de Hipótesis para evaluar si alguna de las dos distribuciones (Weibull y Exponencial) es mejor que otra respecto a su estimación, dado que la mayoría de estos tests no esta diseñados para soportar el rankeo ajustado para datos con censura a la derecha.

Una solución es utilizar el coeficiente de correlacion R^{2} de la regresión lineal de ambas distribución y compararlos con el fin de elegir el mayor, dado que este valor indica cuán explicativo es el modelo de regresión lineal obtenido.

Otra solución es algo más compleja. Luego de obtener los parámetros de la regresión lineal, es posible obtener valores de tiempos ”ficticios” de forma que una estimación de parámetros simple sobre estos tiempos ”ficticios” (que no consideran censura a la derecha) sea la misma que usando los valores de tiempo reales (que si consideran censura a la derecha). De esta forma, los tiempos ”ficticios” podrían contrastarse en un test de hipotesis con la distribución selecciona y obtener los respectivos indicadores acerca de bondad del ajuste.

Por ejemplo, para la distribución Weibull, luego de obtener los parámetros α y β, podemos obtener los tiempos ”ficticios” de esta forma:

ln ( ln (\frac{1}{1 − \frac{r(j)−0,3}{n+0,4}} )) = β ln(t) − β ln(α)

ln(t) = \frac{ln ( ln (_{\frac{1}{1− \frac{r(j)−0,3}{n+0,4}}} )) + β ln(α)}{β}

t = exp[\frac{ln ( ln (_{\frac{1}{1− \frac{r(j)−0,3}{n+0,4}}} )) + β ln(α)}{β}]

en donde reemplazamos r(j) por 1, 2, ..., n donde n es la cantidad arbitraria de tiempos ficticios a elegir. Esta elección consecutiva de valores a reemplazar por r(i) se debe a que el rankeo ajustado perturba los valores j originales (los cuales son números enteros). Finalmente, se obtiene un conjunto n de números que al ajustarlos sin censura a la derecha debería generar una distribución idéntica a la original, y finalmente, usar estos números para el constraste de hipótesis.

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 10

4 APÉNDICE  RMES

4. Apéndice

4.1. Aproximación Ecuación 6

La ecuación es la siguiente

ln ( ln (\frac{1}{1 − \frac{N(t)}{n}} )) = β ln(t) − β ln(α)  (19)

Si el número N(t) es muy pequeño comparado con n, el logaritmo natural interno del lado izquierdo puede aproximarse por el primer término para series de Taylor. Existen muchas funciones con expasión en series conocidas. Dentro ellas, y que no será útiles, están:

ln(1 + x) = ∑_{n=1}^{∞} \frac{(−1)^{n+1}}{n} x^{n} = x − \frac{x^{2}}{2} + \frac{x^{3}}{3} − \frac{x^{4}}{4} + ...  (20)

\frac{1}{1 − x} = ∑_{n=0}^{∞} x^{n} = 1 + x + x^{2} + x^{3} + ...  (21)

Ambas series sólo pueden ser utilizadas si |x| < 1. Si definimos y = N(t)/n, la condición se cumple. Dadas las propiedades lineales de las series de Taylor, es posible resolver ”series sobre series”. Para ello, reescribiremos

ln (\frac{1}{1 − \frac{N(t)}{n}} )

para asemejarla a la ecuacion 20 de la siguiente manera:

ln (\frac{1}{1 − y} )

ln (1 + \overbrace{(\frac{1}{1 − y} − 1)}^{x} )

La ecuación 21 se puede reescribir como:

\frac{1}{1 − y} − 1 = y + y^{2} + y^{3} + ...

Y como sólo deseamos el primer término, la aproximación queda:

\frac{1}{1 − y} − 1 = y  (22)

Ahora, para el primer término de la ecuación 20 obtenemos:

ln(1 + x) = x  (23)

Dado que para la ecuacion 23

x = \frac{1}{1 − y} − 1

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 11 Y según la aproximación hecha en 22 obtenemos que:

ln (1 + (\frac{1}{1 − y} − 1)) = x  (24)

Y remplazando y, finalmente obtenemos:

ln (\frac{1}{1 − \frac{N(t)}{n}} ) ≈ \frac{N(t)}{n}  (25)

4.2. Regresión Lineal Simple

La regresión lineal simple se basa en solucionar el siguiente problema:

mín_{A,B} ∑(y_{i} − ye_{i})^{2}  (26)

donde ye_{i} = Ax_{i} + B, x_{i} es la variable independiente e y_{i} es la variable dependiente.

Lo que se desea es minimizar el error entre el valor de y_{i}, el cual es un dato del problema, respecto a la estimación Ax_{i} − B, donde x_{i} también es un dato del problema.

Teniendo ambos conjuntos de datos \overline{x} e \overline{y} (de igual tamaño) es posible solucionar el problema planteado por la ecuación 26 derivando la función de error parcialmente respecto A y B, igualando a cero y substituyendo el resultado de una ecuación sobre la otra.

Derivando parcialmente para A:

\frac{d}{∂A} ∑(y_{i} − ye_{i})^{2} = \frac{d}{∂A} ∑(y_{i} − Ax_{i} − B)^{2}
= 2 ∑(y_{i} − Ax_{i} − B)\frac{d}{∂A}(y_{i} − Ax_{i} − B)
= 2 ∑(y_{i} − Ax_{i} − B) · −x_{i}

e igualando a cero:

∑(y_{i} − Ax_{i} − B) · −x_{i} = 0
− ∑ x_{i}y_{i} + A ∑ x_{i}^{2} + ∑ Bx_{i} = 0
A = \frac{∑ x_{i}y_{i} − ∑ Bx_{i}}{∑ x_{i}^{2}}

(27)

Derivando parcialmente para B:

\frac{d}{∂B} ∑(y_{i} − ye_{i})^{2} = \frac{d}{∂B} ∑(y_{i} − Ax_{i} − B)^{2}
= 2 ∑(y_{i} − Ax_{i} − B)\frac{d}{∂B} (y_{i} − Ax_{i} − B)
= 2 ∑(y_{i} − Ax_{i} − B) · −1

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 12 e igualando a cero:

∑(y_{i} − Ax_{i} − B) = 0
∑ y_{i} − A ∑ x_{i} − nB = 0
B = \frac{∑ y_{i} − A ∑ x_{i}}{n}

(28)

Uniéndo las ecuaciones obtenidas en 27 y 28, obtenemos que:

A = \frac{n ∑ x_{i}y_{i} − ∑ x_{i} ∑ y_{i}}{n ∑ x_{i}^{2} − ( ∑ x_{i})^{2}}  (29)

B = \frac{∑ y_{i} − A ∑ x_{i}}{n}  (30)

Para el caso de que nuestra ecuación lineal sea aún más simple, o sea, ye_{i} = Ax_{i}, el procedimiento es similar:

Derivando parcialmente para A:

\frac{d}{∂A} ∑(y_{i} − ye_{i})^{2} = \frac{d}{∂A} ∑(y_{i} − Ax_{i})^{2}
= 2 ∑(y_{i} − Ax_{i})\frac{d}{∂A}(y_{i} − Ax_{i})
= 2 ∑(y_{i} − Ax_{i}) · −x_{i}

e igualando a cero:

∑(y_{i} − Ax_{i}) · −x_{i} = 0
− ∑ x_{i}y_{i} + A ∑ x_{i}^{2} = 0
A = \frac{∑ x_{i}y_{i}}{∑ x_{i}^{2}}

(31)

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 13

REFERENCIAS  RMES

Referencias

[1] Weibull Analysis. www.mathpages.com. Documento de libre distribución disponible en http://www.mathpages.com/home/kmath122/kmath122.htm a la fecha 10/01/2011.

[2] Taylor series. Wikipedia. Documento de libre distribución disponible en http://en.wikipedia.org/wiki/Taylor series a la fecha 10/01/2011.

[3] Serie de Taylor. Wikipedia. Documento de libre distribución disponible en http://es.wikipedia.org/wiki/Serie de Taylor a la fecha 10/01/2011.

[4] Distribución Binomial. Wikipedia. Documento de libre distribución disponible en http://es.wikipedia.org/wiki/Distribución binomial a la fecha 11/01/2011.

Centro de Desarrollo de Gestión Empresarial. 1 Oriente 1097, Vina del Mar, Chile. Fono:(56) (32)688987, Fax:(56) (32)2684079 14