Módulo Montecarlo Software MES CGE Consultores

René González Leiva aplikadito@gmail.com

Estudiante Ingeniería Civil Informática Universidad Técnica Federico Santa María Julio, 2008

1. Introducción

La simulación Montecarlo básicamente permite obtener muestreos aleatorios de situaciones reales basado en modelos de comportamiento de las variables que participan en cierto problema.

La implementación hecha en el software MES de la simulación montecarlo permite asignar distintas distribuciones de probabilidad a las variables:

Tiempo entre Fallas (TBF) Tiempo de Reparación (TTR) Costo de Reparación (CR)

Con estas variables, se calculan los siguientes indicadores:

Disponibilidad Costo de Mantención Costo de la Falta

Los resultados son probabilísticos, por lo que se utilizan histogramas para presentar la información.

2. Distribuciones en Módulo Montecarlo para MES

El módulo Montecarlo para MES implementa las siguientes distribuciones de probabilidad: Exponencial Exp(λ) Weibull W eibull(α, β) Normal N(μ, σ^{2}) LogNormal Ln(μ, σ^{2}) Triangular T r(a, b, c) Uniforme U(a, b) Dirac Dirac(p)

Se definen funciones comunes a todas las distribuciones:

Una función que devuelva la función de densidad de la distribución f(x)

Una función que devuelva la función de probabilidad acumulada de la distribución F (x)

Una función que devuelva la función inversa de la función de distribución de probabilidad acumulada F^{−1}(x)

3. Obtención de números aleatorios

Supongamos que se tiene una variable X que se puede modelar con una exponencial de media 0,007, Exp(0, 007). Para obtener valores de la variable X acorde con su comportamiento exponencial, se realiza lo siguiente:

1. Se obtiene un número aleatorio uniformemente distribuido entre 0 y 1 (utilizando la clase Random y el método nextDouble())

2. Luego, con ese número, se llama al método que obtiene la función inversa de la función de distribucion acumulada.

3. El resultado de la inversa es un número dentro del rango definido por la distribución exponencial asociada.

Por ejemplo, supongamos que el número aleatorio obtenido fue 0,4. Sabemos que la función de distribución acumulada de una exponencial es

F (x) = 1 − e^{−(λx)}

La inversa se puede obtener directamente

F (x)^{−1} = −\frac{e^{1−y}}{λ}

donde y es un número entre 0 y 1. De forma genérica

F (x)^{−1} = −\frac{e^{1−U(0,1)}}{λ}

Ahora, simplemente se reemplazando los valores: y = 0,4; λ = 0,007

F (x)^{−1} = 260,3

Con esto, se realiza un mapeo de una función uniforme entre 0 y 1 a una exponencial con media 0.007. Y en general, este procedimiento es similar para toda distribución que posea una función inversa de probabilidad acumulada definida. Este método se conoce como Algoritmo de la Transformada Inversa y se puede definir, de forma general, de la siguiente manera:

Sea U una variable aleatoria uniforme en (0,1). Para cualquier función de distribución continua F, invertible, la variable aleatoria X definida como

X = F^{−1}(U)

tiene distribución F [F^{−1} se define como el valor de x tal que F (x) = u]

4. Obtención de números aleatorios para Normal y LogNormal

4.1. Usando método nextGaussian() de la clase Random

Dado que no es posible aplicar el método anterior a estas distribuciones, se utilizó un generador de números aleatorios normalmente distribuido (utilizando la clase Random y el método nextGaussian()) de manera estandar. Para obtener un número aleatorio respecto a una Normal(μ, σ^{2}), simplemente, se realiza un cambio de variables

n = N(0, 1)
z = σ · n + μ

donde z es un número N(μ, σ)-mente distribuido.

Para obtener un valor aleatorio asociado a una LogNormal con media μ y desviación estandar σ, primero, utilizan las siguientes variables auxiliares

μ_{aux} = lg (\frac{μ^{2}}{\sqrt{μ^{2} + σ^{2}}} )

σ_{aux} = \sqrt{lg (1 + (\frac{σ}{μ})^{2})}

Para después realizar el cambio de variables

z = e^{σ_{aux}·n+μ_{aux}}

donde z es un número Log(μ, σ)-mente distribuido. Este es el método que actualmente se está usando.

4.2. Método de búsqueda recursiva

Este método no se está utilizando actualmente debido a que su performance no es buena. No obstante, su implementanción no ha sido borrada debido a que en el futuro eventualmente podría ser usada. Se adjunta a continuación su explicación

Primero, se define una distribucion Normal Estandar, N(0,1) que sirve de base para los calculos de las distribuciones de probabilidad acumulada y de la inversa. Las razones se detallan a continuación.

Una distribución Normal no posee una función de distribucion acumulada definida. El cálculo para la obtención de esta función se realiza utilizando métodos de integración numérica, o sea, que utilizan la idea básica de Riemann, dividir el dominio de la función en varios intervarlos y calcular el area de estos intervalos con polígonos limitados en su altura por los valores de la función, ya sea por exceso o por defecto. El método se Simpson es uno de los mejores métodos de integracion numérica que hay. Crea polígonos un poco mas complejos a lo largo de la curva y calcula con bastante precisión el area bajo ella. La única dificultad de este método es que necesita de los límites superior e inferior de integración, los cuales no existen para el caso particular de una Normal. Para ello, se utiliza la Normal Estandar, para la cual se conocen limites superiores e inferiores que son factibles de utilizar. En este caso, se utiliza como limite inferior -5 y limite superior 5, dado que las probabilidades de que los números estén fuera del intervalo [-5,5] es ínfima. Por lo tanto, se utilizan estos límites para decir que una Normal Estandar tiene un dominio que sólo existe dentro del intervalo ya definido. Ahora, es posible utilizar el método de Simpson para integrar una Normal Estandar.

Si se desea obtener la función de probabilidad acumulada para una Normal(μ, σ^{2}) para una valor x (es decir, se desea saber cuanto vale la integral desde −∞ hasta el punto x) se puede utilizar el siguiente cambio de variable:

z = \frac{x − μ}{σ^{2}}

donde z es un valor mapeado para el dominio de una Normal Estandar. Luego, se calcula la integral de una Normal Estandar desde −∞ (que en nuestro caso particular hemos definido como -5) hasta el punto z usando el método de Simpson.

De forma idéntica, para una LogNormal se utiliza:

z = \frac{ln(x) − μ}{σ^{2}}

y se calcula la integral para la Normal Estandar usando z.

Finalmente, es posible realizar el cálculo de la función inversa de la distribución acumulada usando un algoritmo recursivo. Para ello, se debe realizar un mapeo inverso. Primero, supongamos que se tiene un número p aleatorio uniformemente distribuido entre 0 y 1. Lo que deseamos saber es con que valor z en la función de distribución acumulada de una Normal Estandar es posible obtener el valor p. En otras palabras, encontrar el valor z tal que la integral entre −∞ y z para una Normal Estandar sea igual a p.

Este valor z se obtiene realizando una búsqueda binaria entre los limites de la Normal Estandar anteriormente definidos de la siguiente manera:

1. Se comienza con limites definidos ( [-5,5] )

2. Se evalua en la mitad de los limites actuales ( F(0) )

3. Si el valor objetivo p es mayor que la evaluación en el punto medio, se corre el limite inferior hasta el punto medio y se vuelve a evaluar en la mitad de los nuevos limites.

4. Si el valor objetivo p es menor que la evaluación en el punto medio, se corre el límite superior hata el punto medio y se vuelve a evaluar en la mitad de los nuevos limites.

5. El algorimo para cuando se la diferencia entre el punto evaluado y la meta p es menor a cierta precisión.

Esto es análogo al proceso de búsqueda binaria de, por ejemplo, un nombre en un directorio telefónico, buscar iterativamente en las mitades de las mitades.

De esta forma, se obtiene un valor z que cumple las condiciones ya mencionadas. Este valor se mapea para el dominio de una Normal(μ, σ^{2}) de la siguiente manera:

x = σ ∗ z + μ

donde x es el valor mapeado para una Normal(μ, σ^{2}) respecto de z en una Normal Estandar. Para una LogNormal(μ, σ^{2}) el cambio de variable es:

x = e^{σ∗z+μ}

donde x es el valor mapeado para una LogNormal(μ, σ^{2}) respecto de z en una Normal Estandar. A pesar de tener definido y resuelto el problema de la obtención de inversas, es necesario definir

cotas inferiores y superiores para las distribuciones, debido al problema computacional de graficar las funciones y además del requerimiento para el módulo de que las distribuciones puedan ser truncadas por cotas definidas por el usuario. Se deben definir cotas superiores e inferiores por defecto. Éstas, arbitrariamente, se definen de la siguiente manera:

Cota Inferior: un valor x inf tal que la probabilidad desde −∞ hasta aquel punto sea 0,1 %. O sea, la funcion inversa de probabilidad acumulada de la distribucion para el valor 0,001.

Cota Superior: el valor x sup tal que la probabilidad desde −∞ hasta aquel punto sea 99,9 %. O sea, la funcion inversa de probabilidad acumulada de la distribucion para el valor 0.999.