Simulacin de Sistemas

1
Simulación de Sistemas

• CO-4311 - Estadística para la Calidad,
  Productividad y Simulación

2
Componentes de un modelo

• Las variables del sistema, las cuales se
  clasifican en variables de entrada
  (independientes o explicativas) y variables de
  salida (dependientes o de respuesta).
• Las relaciones entre las variables, las cuales
  pueden venir dadas por ecuaciones diferenciales
  (funciones de transferencia), por estructuras
  lógicas, etc.

3
Clases de modelos

• Modelos Determinísticos.
• Modelos Estocásticos

Dado que en la mayor parte de los sistemas existe
incertidumbre, los modelos estocásticos son en
general más apropiados. Sin embargo, son también
más complejos y requieren de mayor información
para su implmentación.
4
Por qué simular?

• Calcular las distribuciones de las variables de
  salida puede ser complicado o incluso imposible,
  incluso para distribuciones sencillas.
• Aun cuando sea posible determinar las
  distribuciones, hallar estadísticos relacionados
  con ellas puede ser difícil.

5
Generación de números aleatorios uniformes

• Antiguamente se utilizaban tablas de números
  uniformes.
• Actualmente se utilizan algoritmos
  computacionales.
• Se trata de números pseudoaleatorios, ya que se
  genera una sucesión de números que pasa todas las
  pruebas de aleatoriedad.
• En realidad se genera una distribución discreta
  que aproxima a la verdadera.

6
Generadores Congruneciales Lineales

• Por mucho, los más exitosos generadores de
  números aleatorios que se conocen hoy en día son
  casos particulares de la recursión

De los valores de a, c, m y X0 dependen las
propiedades del generador.
7
Propiedades de los GCL

• La secuencia generada por un GCL es cíclica. El
  número de elementos distintos en la sucesión se
  le denomina período.
• En general se desea que el periodo sea grande.
• El máximo período posible para un GCL es m.
  Puede decir por qué?

8
Propiedades de los GCL (cont)

• Teorema 1.- Un GCL tiene período máximo m si y
  solo si
• c es un primo relativo de m
• b a - 1 es un múltiplo de p para todo primo p
  que sea divisor de m
• b es un múltiplo de 4 si m es un múltiplo de 4.
• Adicionalmente, el valor de a deber ser grande de
  modo que las primeras cifras significativas
  luzcan aleatorias.

9
Propiedades de los GCL (cont)

• Por ejemplo, el generador de Ripley (1987) toma
  la forma
• Verifique que este generador tiene período máximo.

10
Propiedades de los GCL (cont)

• La potencia de un GCL con período máximo esta
  definida como el menor entero s tal que bs mod m
  0.
• Esta propiedad está relacionada con la
  dependencia entre números consecutivos del
  generador, y por lo tanto, es una medida de la
  aletoriedad del mismo.
• En general, se requiere que s gt 4. Cuál es la
  potencia del generador de Ripley?

11
Propiedades de los GCL (cont)

• Para a 85086, c 3, m 425425 (con lo cual s
  3) al graficar 500 pares ordenados de números
  consecutivos se obtiene

12
Generación de variables aleatorias no uniformes

• Todas las técnicas se basan en generadores
  uniformes, y por tanto comparten las
  características de estos.
• El principal método de generación que estudiares
  es el método de inversión. También estudiaremos
  algunas técnicas basadas en transformaciones.

13
El método de inversión

• Teorema 2.- Sea F una distribucíon contínua en R
  con inversa F-1 definida por
• Si UUni(0,1) entonces F-1(U) tiene distribución
  F. Además, si X tiene distribución F, entonces
  F(X) es uniforme en (0,1).

14
El método de inversión (cont)

• La prueba de este teorema se basa en el cálculo
  de una función de una variable aleatoria.
• Añadir el ínfimo y la desigualda en la definición
  de la inversa son necesarios para eliminar
  ambigüedades que se presentan cuando ciertos
  intervalos tienen probabilidad 0 (por ejemplo, en
  distribuciones discretas).

15
El método de inversión (cont)

• Por qué el ínfimo?

Densidad
Distribución
16
El método de inversión (cont)

• Por qué la desigualdad?

Densidad
Distribución
17
Inversión para distribuciones crecientes

• Si la distribución es estrictamente creciente
  entonces no existen ambigüedades en la definición
  de F-1(U), lo que simplifica el algorítmo a
  simplemente evaluar un número aleatorio
  UUni(0,1) en esta función.

18
Inversión para distribuciones crecientes (cont)

• Ejemplo 1.- La distribución exponencial
• La densidad exponencial viene dada por
• y por tanto su distribucón es

19
Inversión para distribuciones crecientes (cont)

• Así pues si UUni(0,1) entonces
• Esto permite generar números que siguen una
  distribución exponencial, si se dispone de número
  que siguen una distribución uniforme en (0,1).

20
Inversión para distribuciones crecientes (cont)
21
Inversión para variables aleatorias discretas

• Supongamos que la variable aleatoria X toma
  valores sobre un conjunto de números consecutivos
  A ? N. En este caso el algoritmo de inversión se
  reduce a
• Generar UUni(0,1).
• Encontrar X tal que F(X - 1) ? U lt F(X).
• El número X así generado sigue la distribución
  dada por F(X).

22
Inversión para variables aleatorias discretas
(cont)

• Ejemplo 2 Distribución de Bernoulli.

23
Inversión para variables aleatorias discretas
(cont)

• Así pues el algorítmo se reduce a
• Generar U Uni(0,1).
• Escoger x 0 si F(-1) ? U lt F(0) ?
• 0 ? U lt 1 - p
• Escoger x 1 si F(0) ? U lt F(1) ?
• 1 - p ? U lt 1
• Hay otras formas adecuadas para escribir estas
  condicones, puede sugerir cuál?

24
Transformaciones

• En muchos casos la variable aleatoria que se
  desea generar pude escribirse como una función de
  otra variable (u otras variables) aleatorias que
  son más fáciles de generar.
• En estos casos, los algorítmos de generación son
  directos.

25
Transformaciones (cont)

• Ejemplo 3.- Distribución Gamma.
• La variable aleatoria X Gamma(n,l) con densidad
  
• no puede generarse directamente a través del
  método de inversión (verifiquelo).

26
Transformaciones (cont)

• Sin embargo, si n es entero, la variable
  aleatoria gamma puede escribirse como una suma de
  n variables aleatorias exponenciales. Esto nos
  sugiere el algoritmo
• Generar n variables aleatorias exponeciales E1,
  E2, .... En con parámetro l.
• Calcular X E1 E2 .... En. Ahora
• X Gamma(n, l).

27
Transformaciones (cont)

• Una versión computacionalmente más eficiente del
  mismo algoritmo es
• Generar U1, U2, .... Un Uni(0,1) e
  independientes.
• Calcular X - l ln(U1U2....Un ). Ahora
• X Gamma(n, l)
• Puede usted decir cómo se obtiene este algoritmo
  en base al anterior y por qué es más eficiente?

28
Transformaciones (cont)

• Ejemplo 4.- Método de Box-Muller (1958) para
  densidades normales.
• Sean X,YN(0,1) y consideremos la transformación
• de modo que la distribución conjunta viene dada
  por

29
Transformaciones (cont)
30
Transformaciones (cont)

• Así, vemos que las variables R y Q son
  independientes con densidades

31
Transformaciones (cont)

• El algoritmo final es
• Generar U1 Uni(0,1) y calcular ? 2p U1.
• Generar U2 Uni(0,1) y calcular
• R (-2ln U2 )1/2.
• Calcular X R cos ?, Y R sen ?. Las variables
  resultantes X e Y son ambas normales estandar.

32
Transformaciones (cont)

• Ejemplo 5.- Distribución uniforme discreta en
  1,k. Se desea generar números enteros entre 1 y
  k de tal modo que todos tengan la misma
  probabilidad. Para ello
• Genere U Uni(0,1).
• Haga X ?kU? 1.
• Entonces X tiene la distribución deseada. Este
  método es más rápido que inversión.

33
Relaciones entre algunas v.a. de uso frecuente.
34
Otras técnicas

• Rechazo. Se generan números aleatorios según una
  distribución similar, y se establece una
  condición para su aceptación.
• Cociente de uniformes. Corresponde a una
  modificación de la técnica de rechazo.
• Monte Carlo Markov Chains (MCMC).

35
Distribuciones truncadas

• Son aquellas en las que el soporte natural de la
  distribución se reduce. Por ejemplo, una
  Normal(0,1) que se restringe al intervalo (-2,2).
  La densidad resultante es proporcional a la
  original.
• Para simularla, se generan números según la
  densidad original y se descartan aquellos que
  caigan fuera del soporte truncado (rechazo).

36
Variable aleatorias mixtas

• Son aquellas que presentan puntos de masa en
  algunos valores y en el resto corresponden a una
  v.a. continua.
• La generación se hace en dos etapas primero se
  decide si estamos en un punto de masa o en una
  región continua, y luego se genera localmente
  dependiendo de lo anterior.

37
Variable aleatorias mixtas (cont)

• Ejemplo 6.- Un cierto componente de repuesto
  tiene una probabilidad p 0,05 se venir
  defectuoso del proveedor. Cuando el componente
  no tiene problemas de manufactura, el tiempo
  durante el cual funciona bien sigue una densidad
  exponencial con media 1 año.
• Puede dibujar la densidad correspondiente?

38
Variable aleatorias mixtas (cont)

• El tiempo de vida de este componente es un
  ejemplo de una variable aleatoria mixta con un
  punto de masa en T 0. El algoritmo para
  simularlo es el siguiente
• Generar U1 Uni(0,1).
• Si U1 ? 0.05 entonces T 0. Si U1 gt 0.05
  entonces generar U2 Uni(0,1) y hacer
• T - ln (1 - U2).
• Puede plantear un algoritmo alternativo?

39
Generación de distribuciones empíricas

• A veces los datos correspondientes a alguna de
  las variables no se ajustan adecuadamente a
  ninguna distribución conocida. Algunas
  alternativas son
• Remuestreo (bootstrap) si la muestra es grande.
• Generación a partir del histograma.
• Estimación de densidades (por ejemplo, a través
  de combinaciones lineales de kernels).

40
Generación de distribuciones empíricas

• En el remuestreo se escoge un subconjunto de los
  datos originales.
• Para generar a partir del histograma
• Se decide primero a cual categoría pertenece el
  dato, generando de una distribución discreta
  donde la probabilidad de cada categoría es igual
  a su frecuencia observada.
• Se genera X Uni(LI,LS), donde LI y LS son los
  valores que delimitan la categoría.

41
Aplicaciones

• Ejemplo 7.- Cierto equipo utiliza tres
  componentes conectados en serie. El tiempo de
  vida (en meses) para cada uno de ellos viene dado
  por las siguientes distribuciones
• T1 Exponencial (25)
• T2 Weibull (320)
• T3 Lognormal(2,50,6)
• Lo cual quiere decir que las densidades
  correspondientes son

42
Aplicaciones (cont)

• Puede decir como generar números aleatorios que
  sigan cada una de estas distribuciones?

43
Aplicaciones (cont)

• Al estar los componentes conectados en serie el
  tiempo de vida del sistema Ts es
• Puede ver el por qué de esta relación?
• Al simular el comportamiento de 10.000 sistemas
  de este tipo se obtiene

44
Aplicaciones (cont)

• Histograma y diagrama de caja del tiempo de vida
  del sistema.
• Q1 5.158 Q2 8.301 Q3 12.210 ET 9.000

45
Aplicaciones (cont)

• La distribución es asimétrica.
• El tiempo medio de vida del sistema completo es
  de 9 horas.
• A pesar de que en algunos casos el sistema puede
  durar incluso más de 20 horas, estos componentes
  son muy raros, ya que menos del 25 excede las 13
  horas de vida.

46
Aplicaciones (cont)

• Probabilidad de falla en función del tiempo.
• Puede dar un tiempo de garantía para el sistema?

47
Aplicaciones (cont)

• Proporción de fallas debidas a cada componente.

48
Aplicaciones (cont)

• Claramente, el componente 1 es la segunda causa
  más importante de fallas, no lejos del componente
  2 . Sin embargo
• ET1 25,02 ET2 17,75 ET3 14,53
• Es decir, es por mucho el componente con el
  tiempo de vida promedio más largo, Puede
  explicar esta paradoja?

49
Aplicaciones (cont)

• Simular un sistema con tres componentes idénticos
  y verificar que los tres producen la falla con la
  misma frecuencia.
• También sería interesante determinar si alguno de
  los componentes causa las fallas más tempranas
  y si alguno las más tardías.

50
Aplicaciones (cont)

• Ejemplo 8.- Una tarjeta de circuito impreso
  tiene un cierto número de huecos que se hacen
  utilizando un taladro numérico controlado por una
  computadora. La computadora tiene un número de
  fallas por tarjeta aleatorio con distribución
  Poisson (6), y si la computadora falla, la
  probabilidad de que el taladro no haga hueco es
  0,15.

51
Aplicaciones (cont)

• Sea F el número de fallos por tarjeta que comete
  la computadora, y H el número de huecos no hechos
  por el taladro en una tarjeta. Entonces

52
Aplicaciones (cont)
Simulando a partir de estas distribuciones se
obtienen las siguientes marginales
53
Aplicaciones (cont)

• El número de huecos faltantes por tarjeta es
  menor al número de fallos de la computadora.
• La probabilidad de tener una tarjeta operativa
  (sin huecos faltantes) es apenas de 0.399, lo que
  significa una proporción de defectos altísima.
  Puede sugerir una forma de corregir el problema?

54
Aplicaciones (cont)

• Ejemplo 9.- Una planta que trabaja solo bajo
  pedidos desea determinar si puede satisfacer uno
  de 17 unidades para el próximo día habil sin usar
  sobretiempo. La planta trabaja un número de
  horas/día T que sigue una distribución uniforme
  en (78), y el número de piezas que produce sigue
  una distribución de Poisson con parámetro 2,5T
  piezas/dia.

55
Aplicaciones (cont)

• Una forma de responder podría ser calcular el
  número promedio de piezas que se pueden fabricar
  en un día.
• En apariencia podemos cumplir con el pedido. Es
  eso verdad?

56
Aplicaciones (cont)

• Del enunciado es claro que
• Simulando el comportamiento del sistema durante
  10.000 días podemos obtener la distribución

57
Aplicaciones (cont)

• que no es más que la probabilidad de no fabricar
  en un día suficientes piezas para un pedido de
  tamaño k.

58
Aplicaciones (cont)

• Es claro que al menos un 31 de las veces no
  podremos cumplir con el pedido.
• Si el tamaño del pedido aumenta a 19 unidades
  (solo ligeramente sobre el promedio), la
  probabilidad es del 49.
• El problema podría haber sido incluso más grave
  si la distribución hubiese tenido sesgo
  contrario.

59
Aplicaciones (cont)

• Este ejemplo muestra la diferencia entre una
  planta de producción continua y una bajo pedido
  (JIT, teoría de restricciones).
• Normalmente la estimación del número de piezas a
  fabricar en un día proviene a su vez de otra
  simulación.
• Que hubiese pasado con un pedido más grande,
  digamos de 150 unidades en 8 días hábiles? Por
  qué?.

60
Tamaño de la Simulación

• Un aspecto fundamental es cuántas iteraciones
  necesitamos en nuestra simulación. Está claro
  que, en términos generales, entre más mejor. Sin
  embargo, existe un tope a partir del cual las
  mejoras que se obtienen en la estimación no
  compensan el esfuerzo realizado.

61
Tamaño de la Simulación (cont)

• En muchos casos la experiencia práctica, el
  sentido común y/o la disponibilidad de tiempo y
  recursos nos pueden dar una cotas para n.
• Tiempo durante el cual operará la planta.
• Número de artículos de una determinada serie que
  se espera vender antes de modificarla.
• Tiempo que tarda correr la simulación.

62
Tamaño de la Simulación (cont)

• Teorema 3.- Ley de los Grandes Números para
  proporciones. Sea p la probabilidad del evento
  A, y fA nA/n la frecuencia observada de A.
  Entonces, para cualquier número positivo ? se
  tiene

63
Tamaño de la Simulación (cont)

• Escoger ?, la discrepancia máxima que vamos a
  permitir, y ?, la cota superior de la
  probabilidad de la discrepancia ?.
• Si no se conoce p (lo más común), el número de
  simulaciones n viene dado por

64
Tamaño de la Simulación (cont)

• Ejemplo 10.- Recuerde el ejemplo 6, donde se
  quería obtener la proporción de tarjetas
  defectuosas. Supongamos ? 0,01 (como máximo un
  error en la segunda cifra significativa) y que ?
  0,05 (se desea una probabilidad baja de estar
  lejos de la verdadera probabilidad)

65
Tamaño de la Simulación (cont)

• Teorema 4.- Ley de los Grandes Números para
  promedios. Sean X1,...,Xn variables aleatorias
  iid con E(Xi)? y Var(Xi) ?2. Entonces para
  para cualquier número positivo ? se tiene

66
Tamaño de la Simulación (cont)

• Escoger ?, la discrepancia máxima que vamos a
  permitir, y ?, la cota superior de la
  probabilidad de la discrepancia ?. Entonces
• Como en general ?2 no se conoce por lo que se
  utiliza un estimador. Esto sugiere un algoritmo
  iterativo, ya que la exactitud del estimador
  también depende de n.

67
Simulaciones Dinámicas

• Hasta ahora hemos realizado simulaciones
  estáticas (donde las variables aleatorias no
  dependen del tiempo).
• Para simulaciones dinámicas es importante definir
  dos conceptos
• Estados características relevantes del sistema
  (ej. el número de individuos en una cola).
• Eventos acciones que pueden producir cambios de
  estado (ej. llegada de una persona a la cola).

68
Simulaciones Dinámicas

• Cómo simular procesos dinámicos? Una idea
  atrayente es discretizar el tiempo, y luego
  avanzar en él verificando si en cada instante se
  produce o no un evento.
• Esta estrategia tiene algunos problemas la
  resolución de la discretización es arbitraria,
  los tiempos de ejecución largos, etc.

Tiempo
69
Simulación de Eventos Discretos

• En el caso de sistemas con eventos discretos
  (aquellos cuyos eventos pueden suceder solo en
  tiempos discretos) el truco es llevar una lista
  con los tiempos en los cuales sucederán los
  próximos eventos (LTE). Esto permite realizar
  grandes saltos en el tiempo.

70
Simulación de eventos discretos (cont)

• Para sistemas que cambian de estado
  continuamente, como los sistemas químicos o
  aquellos dirigidos por ecuaciones diferenciales,
  esta metodología no puede ser utilizada.
• Además de la lista de tiempo de eventos futuros
  es necesario rastrear el estado del sistema (ES).

71
Simulación de eventos discretos (cont)
72
Simulación de eventos discretos (cont)

• Ejemplo 9.- Línea de espera con un solo
  servidor. Se desea simular una línea de espera
  en la cual los tiempos en minutos entre llegadas
  de nuevos individuos A Exp(4), mientras que los
  tiempos de servicio en minutos D Exp(3).
  Suponga que el sistema funciona durante 8 horas,
  que no hay límites para el tamaño de la cola y
  que al inicio no hay nadie en la fila.

73
Simulación de eventos discretos (cont)

• En este caso el estado del sistema está dado por
  el número de personas n que se encuentran dentro
  de él (cola servicio).
• En cuanto a los eventos, en este caso son de solo
  dos tipos llegadas y salidas. Así LE
  (tA,tB).
• Dependiendo del objetivo concreto también es
  necesario mantener algunos contadores o
  acumuladores.

74
Simulación de eventos discretos (cont)

• Asignar t 0
• Asignar n 0
• Generar tA Exp(4)
• Asignar tD ?
• Repetir mientras min(tA, tD) lt 480
• Si tA lt tD
• Asignar t tA
• Asignar n n 1
• Generar ?tA Exp(4)
• Asignar tA t ?tA
• Si n 1
• Generar ?tD Exp(3)
• Asignar tD t ?tD

• En caso contrario
• Asignar t tD
• Asignar n n - 1
• Si n 0
• Asignar tD ?
• En caso contrario
• Generar ?tD Exp(3)
• Asignar tD t ?tD

75
Simulación de eventos discretos (cont)

• Si los sujetos en cola al final del día quedan en
  la misma hasta que se reinicien las actividades
  (e.j. productos en una fábrica), este código es
  suficiente.
• En caso contrario hay que añadir
• Repetir mientras n gt 0
• Asignar t tD
• Asignar n n - 1
• Generar ?tD Exp(3)
• Asignar tD t ?tD

76
Simulación de eventos discretos (cont)

• Si se tratase de una cola de atención al público,
  cómo fijaría usted la jordana de trabajo diaria
  para sus empleados y los horarios de atención?
• Si se tratase de una fábrica, cómo calcularía la
  capacidad y los requisitos de materia prima?
  Considere todos los casos que crea convenientes.

77
Simulación de eventos discretos (cont)

• Ejemplo 10 Vamos a considerar el caso de una
  taquilla de atención al público a la cual llegan
  en promedio 2 personas cada minuto, mientras que
  el servidor atiende en promedio tres personas
  cada minuto. Si los tiempos entre llegadas y
  entre servicios siguen una distribución
  exponencial cómo podría determinar la hora de
  cierre óptima?

78
Simulación de eventos discretos (cont)

• El funcionamiento de este sistema es
  esencialmente el mismo que el del ejemplo 9,
  ahora con A Exp(1/2), mientras que D Exp(1/3)
  (recuerde que estamos usando a la esperanza como
  parámetro de la distribución exponencial).
• Ahora bien, qué significa una hora de cierre
  óptima?

79
Simulación de eventos discretos (cont)

• El primer costo que en el que podemos incurrir es
  un costo de sobretiempo por ley, la jornada
  laboral es de 8 horas diarias, y en caso de
  excederse, se han de cancelar horas extras. Este
  costo es fácil de determinar en base al sueldo
  mensual del cajero y a la ley del trabajo.
  Puede decir cuál es la solución óptima si solo
  se toma en cuenta este costo?

80
Simulación de eventos discretos (cont)

• El problema se debe a que la función de costos
  que hemos definido hasta ahora es monótona y no
  cóncava.
• Para resolver esto necesitamos adicionar un
  término de costos que tienda a decrecer según se
  incremente el tiempo de la jornada laboral. Este
  término está relacionado con la satisfacción del
  cliente, la cual es difícil de medir
  monetariamente

81
Simulación de eventos discretos (cont)

• Para saltarnos este problema es entonces
  razonable estimar una distribución para el costo
  de sobretiempo en cada jornada posible y escoger
  en base a nuestro conocimiento del problema
  aquella donde el costo promedio por horas extra
  tome un valor razonable.

82
Simulación de eventos discretos (cont)

• colaunservidor_function(r,secuencia)
• cu_200
• costopromedio_rep(0,length(secuencia))
• for(j in 1length(secuencia))
• costo_rep(0,r)
• for(i in 1r)
• t_0
• n_0
• ta_rexp(1,1/4)
• td_ta1000
• while(min(ta,td)ltsecuenciaj)
• if(talttd)
• t_ta
• n_n1
• dta_rexp(1,1/4)
• ta_tdta
• if(n1)
• dtd_rexp(1,1/3)

• td_tdtd
• else
• t_td
• n_n-1
• if(n0)
• td_ta1000
• else
• dtd_rexp(1,1/3)
• td_tdtd
• while(ngt0)
• t_td
• n_n-1
• dtd_rexp(1,1/3)
• td_tdtd
• costoi_ifelse(tgt480,cu(t-480),0)
• costopromedioj_mean(costo)
• return(costopromedio)

83
Simulación de eventos discretos (cont)

• Usando el anterior código de R se puede generar
  una curva de costo promedio, la cual luce así.

84
Simulación de eventos discretos (cont)

• Si bien una función de costo es la mejor opción,
  es muy común en ingeniería que la misma sea muy
  difícil de obtener, o que la misma sea tan
  subjetiva como las elecciones que hemos hecho
  hasta ahora para determinar óptimos.

85
Simulación de eventos discretos (cont)

• Modelo gráfico de una cola con un servidor usando
  Extend.

86
Simulación de eventos discretos (cont)