UNIVERSIDADE DA CORUÑA
FACULTAD DE INFORMÁTICA Departamento de Matemáticas
PROYECTO DE FIN DE CARRERA DE INGENIERÍA INFORMÁTICA
Diseño e implementación de una herramienta para la enseñanza y el aprendizaje de la teoría de colas
Autor :
Jorge L. Vega Valle
Director : Ricardo Cao Abad
A Coruña, Enero de 2004
Título del proyecto:
Diseño e implementación de una herramienta para la enseñanza y aprendizaje de la teoría de colas.
Clase de proyecto:
Clásico de ingeniería
Nombre del autor:
Jorge L. Vega Valle
Nombre del director: Ricardo Cao Abad
Miembros del Tribunal
Fecha de lectura y defensa Calificación obtenida
A mi madre
AGRADECIMIENTOS Se me hace bastante complicado el intentar condensar en una sola página los agradecimientos a todas aquellas personas que han tenido una influencia directa o indirecta en la realización de este proyecto. En primer lugar, a mi familia, en especial a mi abuela, mi madre y mi hermana, por su preocupación, su apoyo en los momentos difíciles y la confianza depositada en mí. También me gustaría agradecer a todos aquellos profesores de la facultad, especialmente a mi director de proyecto, Ricardo Cao, por su dedicación, su disposición y por toda la ayuda prestada y a José Mª Domínguez. No sería justo que me quedase sin agradecer tanto a todos los compañeros de prácticas que he tenido a lo largo de estos años: Pablo, Alberto, Víctor, Miguel,... y sobre todo a Juanda, como a todos los compañeros y amigos que he conocido en la facultad. Por último, aunque no por ello menos importante, me gustaría también dar las gracias a mis amigos de Santa María del Mar, especialmente a Juan y a Luis, por su preocupación y las veces que nos hemos reído todos juntos.
RESUMEN Mediante la realización de este proyecto se pretende desarrollar una aplicación que permita resolver los distintos problemas de teoría de colas (modelos y redes de colas) que se puedan plantear, y que al mismo tiempo pueda ser empleada en un ámbito docente. Para ello se ha desarrollado en MATLAB ® una aplicación que resuelve dichos problemas desde una perspectiva dual: aquellos modelos en los cuales tanto la distribución del tiempo entre llegadas como la distribución del tiempo de servicios sean de carácter exponencial, se resolverán de forma analítica, proporcionando una respuesta exacta a dicho problema implementado las distintas fórmulas existentes; por otra parte cuando la distribución del tiempo entre llegadas o la del tiempo de servicios sean o no exponenciales, se resolverá dicho problema empleando técnicas de simulación. En ambos casos el usuario tiene la posibilidad de conocer las características del problema introducido mediante la visualización (gráfica o numérica) de un determinado número de parámetros como los siguientes: número medio de clientes en el sistema, número medio de clientes en la cola, tiempo medio que un cliente está en el sistema, tiempo medio de espera en la cola, eficiencia del sistema, intensidad de tráfico del sistema, probabilidad de que haya un determinado número de clientes en el sistema, probabilidad de que haya un determinado número de clientes en el sistema justo cuando una nueva llegada se está produciendo o la probabilidad de que un determinado cliente esté en la cola o en el sistema durante un tiempo menor que uno establecido. Con el desarrollo de esta aplicación, se pretende ofrecer al usuario la posibilidad de conocer mejor las características del problema introducido e incluso tomar decisiones sobre el mismo si se trata de un problema real.
PALABRAS CLAVE Fase de estabilización, Fase de simulación, Modelo M/M/-, Redes de Jackson abiertas, Redes de Jackson cerradas, Simulación estocástica, Sistemas de espera, Tasa de llegada, Tasa de servicio, Teoría de colas.
Lista de Figuras Figura 1- Modelo de una cola simple .........................................................................15 Figura 2- Elementos existentes en un modelo de colas..............................................19 Figura 3- Red de Jackson abierta con K = 3 ...............................................................28 Figura 4- Red de Jackson cerrada con K = 3 ..............................................................30 Figura 5- Resumen de tipos de simulación ................................................................36 Figura 6- Comportamiento de la simulación..............................................................38 Figura 7- Ciclo de vida para la parte de resolución analítica .....................................42 Figura 8- Ciclo de vida para la parte de simulación...................................................43 Figura 9- Comportamiento gráfico de la simulación..................................................68 Figura 10- Actualización de los valores de "c" y "d" en las redes cerradas ................82 Figura 11- Localización de la carpeta en donde se encuentra la aplicación.............125 Figura 12- Ventana inicial de la aplicación..............................................................126
Índice Capítulo 1: Introducción........................................................................................... 11
1.1- Características iniciales del proyecto ..........................................................12
1.1.1- Justificación del proyecto.........................................................................12 1.1.2- Objetivos del proyecto .............................................................................12 1.1.3- Herramienta utilizada para el desarrollo ..................................................13 1.2- Introducción a la Teoría de Colas ...............................................................14
1.2.1- El srcen...................................................................................................14 1.2.2- Definiciones iniciales...............................................................................15 1.2.3- Objetivos de la Teoría de Colas ...............................................................16 1.2.4- Elementos existentes en un modelo de colas ...........................................17 1.2.5- Terminología y notación ..........................................................................20 1.2.5.1- Terminología.........................................................................................20 1.2.5.2- Conceptos básicos .................................................................................21 1.2.5.3- Notación de Kendall..............................................................................22 1.2.6- Redes de colas..........................................................................................24 1.2.6.1- Introducción a las redes de colas...........................................................24 1.2.6.2- Redes de Jackson abiertas .....................................................................26 1.2.6.3- Redes de Jackson cerradas ....................................................................28 1.3- Introducción a la Simulación ....................................................................... 31
1.3.1- Conceptos básicos....................................................................................31 1.3.2- Experimentación real y simulación..........................................................32 1.3.3- Tipos de simulación .................................................................................33 1.3.3.1- Simulación estática y dinámica.............................................................33 1.3.3.2- Simulación por eventos y por cuantos ..................................................34 1.3.4- Problemas de estabilización y dependencia .............................................36 1.3.5Comportamiento de la simulación 1.3.6- Aplicaciones comunes y uso de la ...........................................................37 simulación.........................................38
Capítulo 2: Aspectos de implementación ..................................................................40
2.1- Características generales de la aplicación .................................................. 41
2.1.1- Filosofía de la implementación ................................................................41 2.1.2- Ciclo de vida empleado............................................................................42 2.1.3de entrada y salida.................................................................43 2.1.4- Parámetros Distribuciones de probabilidad elegidas ..................................................48 2.2- Resolución analítica ......................................................................................50
2.2.1- Características generales ..........................................................................50 2.2.2- Resolución de los modelos.......................................................................50 2.2.2.1- Modelo M/M/1 ......................................................................................50 2.2.2.2- Modelo M/M/s .......................................................................................51 2.2.2.3- Modelo M/M/1/K...................................................................................52 2.2.2.4- Modelo M/M/s/K ...................................................................................54 2.2.2.5- Modelo M/M/1/∞ /H .............................................................................55 2.2.2.6- Modelo M/M/s/∞ /H ..............................................................................57 2.2.2.7- Modelo M/M/s/∞ /H con Y repuestos....................................................58 2.2.2.8- Modelo M/M/∞......................................................................................59 2.2.2.9- Redes de Jackson abiertas .....................................................................60 2.2.2.10- Redes de Jackson cerradas ..................................................................63 2.3- Resolución por simulación............................................................................67
2.3.1- Características generales ..........................................................................67 2.3.2- Implementación de los distintos modelos ...............................................68 2.3.2.1- Bucle de simulación del modelo G/G/1 ................................................69 2.3.2.2- Bucle de simulación del modelo G/G/s.................................................70 2.3.2.3- Bucle de simulación del modelo G/G/1/K ............................................72 2.3.2.4- Bucle de simulación del modelo G/G/s/K .............................................73 2.3.2.5- Bucle de simulación del modelo G/G/1/∞ /H .......................................74 2.3.2.6- Bucle de simulación del modelo G/G/s/∞ /H........................................75 2.3.2.7- Bucle de simulación del modelo G/G/s/∞ /H con Y repuestos..............77 2.3.2.8- Bucle de simulación del modelo G/G/∞ ...............................................78 2.3.2.9 Bucle de simulación de las redes cerradas .............................................80 2.3.2.10- Bucle de simulación de las redes abiertas...........................................83
2.4- Detalles de la implementación ...................................................................... 87
2.4.1- Cálculo de W(t) en los modelos en los que no hay una solución analítica ............................................................................................................................87 2.4.2- Problemas de desbordamiento en el cálculo ............................................88 2.4.3- Problemas de precisión en el cálculo .......................................................89
Capítulo 3: Evaluación de la aplicación................................................................... 91
3.1- Validación de los resultados de la simulación ............................................92 3.1.1- Variación del parámetro de estabilización y del número de clientes de la simulación ..........................................................................................................92 3.1.2- Variación de la intensidad de tráfico........................................................95 3.1.3- Variación del modelo de colas .................................................................97 3.2- Validación de la simulación de las distribuciones de probabilidad........100 3.3- Influencia de la variación de algunos parámetros en el tiempo de ejecución de la simulación .................................................................................111 3.3.1- Variación de la distribución de probabilidad .........................................111 3.3.2- Variación de la intensidad de tráfico......................................................113 3.3.3- Variación del parámetro de estabilización y número de clientes de la simulación ........................................................................................................114 3.3.4- Variación del modelo de colas ...............................................................116
Capítulo 4: Conclusiones y trabajo futuro .............................................................118
Apéndice A: Características de las de distribuciones usadas en la aplicación . 121 Apéndice B: Instalación y ejecución de la aplicación..........................................124 Apéndice C: Contenido del CD ............................................................................. 127
Bibliografía..............................................................................................................128
Capítulo 1: Introducción “No importa en qué cola se sitúe: La otra siempre avanzará más rápido” (Primera Ley de Harper)
“Y si se cambia de cola, aquélla en la que estaba al principio empezará a ir más deprisa” (Segunda Ley de Harper)
Capítulo 1: Introducción
Jorge L. Vega Valle
1.1- Características iniciales del proyecto 1.1.1- Justificación del proyecto Las “colas” son un aspecto de la vida moderna que nos encontramos continuamente en nuestras actividades diarias. En el contador de un supermercado, accediendo a Internet,... el fenómeno de las colas surge cuando unos recursos compartidos necesitan ser accedidos para dar servicio a un elevado número de trabajos o clientes. El estudio de las colas es importante porque proporciona tanto una base teórica del tipo de servicio que podemos esperar de un determinado recurso, como la forma en la cual dicho recurso puede ser diseñado para proporcionar un determinado grado de servicio a sus clientes. Debido a lo comentado anteriormente, se plantea como algo muy útil el desarrollo de una herramienta que sea capaz de dar una respuesta sobre las características que tiene un determinado modelo de colas.
1.1.2- Objetivos del proyecto Estos son algunos de los objetivos que se pretenden alcanzar con el desarrollo de la aplicación: Facilitar al usuario una herramienta que le permita obtener de forma sencilla (tanto gráfica como analíticamente), modelo de colas.
12
las características de un determinado
Capítulo 1: Introducción
Jorge L. Vega Valle
Fomentar el estudio y aprendizaje de la teoría de colas en un ámbito docente, usando la herramienta como un posible complemento de las clases teóricas. Poder estudiar como influye la variación de los parámetros de un modelo sobre el modelo inicial para así, facilitar la toma de decisiones.
1.1.3- Herramienta utilizada para el desarrollo La programación de los distintos modelos de colas conlleva la codificación de un gran número de fórmulas matemáticas, no todas ellas sencillas, que dan respuesta a las características del mismo. Asimismo, el cálculo vectorial y matricial está presente en una gran parte de estas fórmulas matemáticas (por ejemplo: las probabilidades para ir de un nodo a otro forman una matriz, la evolución de (t) a lo largo del tiempo se puede guardar en un vector,...). Por todo ello, y también debido a su facilidad de uso, a la productividad respecto a otros entornos de desarrollo y a que se trata de una herramienta ampliamente extendida en el ámbito docente que puede ayudar a fomentar la enseñanza y el aprendizaje de la teoría de colas se ha elegido MATLAB ® como la herramienta más adecuada para resolver un problema de estas características.
13
Capítulo 1: Introducción
Jorge L. Vega Valle
1.2- Introducción a la Teoría de Colas En muchas ocasiones en la vida real, un fenómeno muy común es la formación de colas o líneas de espera. Esto suele ocurrir cuando la demanda real de un servicio es superior a la capacidad que existe para dar dicho servicio. Ejemplos reales de esa situación son: los cruces de dos vías de circulación, los semáforos, el peaje de una autopista, los cajeros automáticos, la atención a clientes en un establecimiento comercial, la avería de electrodomésticos u otro tipo de aparatos que deben ser reparados por un servicio técnico, etc. Todavía más frecuentes, si cabe, son las situaciones de espera en el contexto de la informática, las telecomunicaciones y, en general, las nuevas tecnologías. Así, por ejemplo, los procesos enviados a un servidor para ejecución forman colas de espera mientras no son atendidos, la información solicitada, a través de Internet, a un servidor Web puede recibirse con demora debido a congestión en la red o en el servidor propiamente dicho, podemos recibir la señal de líneas ocupadas si la central de la que depende nuestro teléfono móvil está colapsada en ese momento, etc.
1.2.1- El srcen El srcen de la Teoría de Colas está en el esfuerzo de Agner Kraup Erlang (Dinamarca, 1878 - 1929) en 1909 para analizar la congestión de tráfico telefónico con el objetivo de cumplir la demanda incierta de servicios en el sistema telefónico de Copenhague. Sus investigaciones acabaron en una nueva teoría denominada teoría de colas o de líneas de espera. Esta teoría es ahora una herramienta de valor en negocios debido a que un gran número de problemas pueden caracterizarse, como problemas de congestión llegada-salida.
14
Capítulo 1: Introducción
Jorge L. Vega Valle
1.2.2- Definiciones iniciales La teoría de colas es el estudio matemático del comportamiento de líneas de espera. Esta se presenta, cuando los “clientes” llegan a un “lugar” demandando un servicio a un “servidor”, el cual tiene una cierta capacidad de atención. Si el servidor no está disponible inmediatamente y el cliente decide esperar, entonces se forma la línea de espera. Una cola es una línea de espera y la teoría de colas es una colección de modelos matemáticos que describen sistemas de línea de espera particulares o sistemas de colas. Los modelos sirven para encontrar un buen compromiso entre costes del sistema y los tiempos promedio de la línea de espera para un sistema dado. Los sistemas de colas son modelos de sistemas que proporcionan servicio. Como modelo, pueden representar cualquier sistema en donde los trabajos o clientes llegan buscando un servicio de algún tipo y salen después de que dicho servicio haya sido atendido. Podemos modelar los sistemas de este tipo tanto como colas sencillas o como un sistema de colas interconectadas formando una red de colas. En la siguiente figura podemos ver un ejemplo de modelo de colas sencillo. Este modelo puede usarse para representar una situación típica en la cual los clientes llegan, esperan si los servidores están ocupados, son servidos por un servidor disponible y se marchan cuando se obtiene el servicio requerido.
Servidor(es) Llegadas Salidas Posiciones de espera Figura 1- Modelo de una cola simple
15
Capítulo 1: Introducción
Jorge L. Vega Valle
El problema es determinar qué capacidad o tasa de servicio proporciona el balance correcto. Esto no es sencillo, ya que un cliente no llega a un horario fijo, es decir, no se sabe con exactitud en que momento llegarán los clientes. También el tiempo de servicio no tiene un horario fijo. Los problemas de “colas” se presentan permanentemente en la vida diaria: un estudio en EEUU concluyó que, por término medio, un ciudadano medio pasa cinco años de su vida esperando en distintas colas, y de ellos casi seis meses parado en los semáforos.
1.2.3- Objetivos de la Teoría de Colas Los objetivos de la teoría de colas consisten en: Identificar el nivel óptimo de capacidad del sistema que minimiza el coste global del mismo. Evaluar el impacto que las posibles alternativas de modificación de la capacidad del sistema tendrían en el coste total del mismo. Establecer un balance equilibrado (“óptimo”) entre las consideraciones cuantitativas de costes y las cualitativas de servicio. Hay que prestar atención al tiempo de permanencia en el sistema o en la cola: la “paciencia” de los clientes depende del tipo de servicio específico considerado y eso puede hacer que un cliente “abandone” el sistema.
16
Capítulo 1: Introducción
Jorge L. Vega Valle
1.2.4- Elementos existentes en un modelo de colas Fuente de entrada o población potencial: Es un conjunto de individuos
(no necesariamente seres vivos) que pueden llegar a solicitar el servicio en cuestión. Podemos considerarla finita o infinita. Aunque el caso de infinitud no es realista, sí permite (por extraño que parezca) resolver de forma más sencilla muchas situaciones en las que, en realidad, la población es finita pero muy grande. Dicha suposición de infinitud no resulta restrictiva cuando, aún siendo finita la población potencial, su número de elementos es tan grande que el número de individuos que ya están solicitando el citado servicio prácticamente no afecta a la frecuencia con la que la población potencial genera nuevas peticiones de servicio. Cliente: Es todo individuo de la población potencial que solicita servicio.
Suponiendo que los tiempos de llegada de clientes consecutivos son 0< t1
Capacidad de la cola: Es el máximo número de clientes que pueden estar haciendo cola (antes de comenzar a ser servidos). De nuevo, puede suponerse
finita o infinita. Lo más sencillo, a efectos de simplicidad en los cálculos, es suponerla infinita. Aunque es obvio que en la mayor parte de los casos reales la
17
Capítulo 1: Introducción
Jorge L. Vega Valle
capacidad de la cola es finita, no es una gran restricción el suponerla infinita si es extremadamente improbable que no puedan entrar clientes a la cola por haberse llegado a ese número límite en la misma. Disciplina de la cola: Es el modo en el que los clientes son seleccionados
para ser servidos. Las disciplinas más habituales son: La disciplina FIFO (first in first out), también llamada FCFS (first come first served): según la cual se atiende primero al cliente que antes haya llegado. La disciplina LIFO (last in first out), también conocida como LCFS (last come first served) o pila: que consiste en atender primero al cliente que ha llegado el último.
La RSS (random selection of service), o SIRO (service in random order), que selecciona a los clientes de forma aleatoria. La disciplina RR (round robin), según la cual se otorga un pequeño cuanto de tiempo de servicio a cada cliente de forma secuencial. Esto viene a equivaler a repartir los recursos de forma igualitaria entre todos los clientes en espera y, por supuesto sólo tiene sentido en algunas circunstancias (como el ámbito de la informática). Mecanismo de servicio: Es el procedimiento por el cual se da servicio a
los clientes que lo solicitan. Para determinar totalmente el mecanismo de servicio debemos conocer el número de servidores de dicho mecanismo (si dicho número fuese aleatorio, la distribución de probabilidad del mismo) y la distribución de probabilidad del tiempo que le lleva a cada servidor dar un servicio. En caso de
18
Capítulo 1: Introducción
Jorge L. Vega Valle
que los servidores tengan distinta destreza para dar el servicio, se debe especificar la distribución del tiempo de servicio para cada uno. La cola, propiamente dicha, es el conjunto de clientes que hacen espera, es
decir los clientes que ya han solicitado el servicio pero que aún no han pasado al mecanismo de servicio. El sistema de la cola : es el conjunto formado por la cola y el mecanismo
de servicio, junto con la disciplina de la cola, que es lo que nos indica el criterio de qué cliente de la cola elegir para pasar al mecanismo de servicio. Estos elementos pueden verse más claramente en la siguiente figura:
Disciplina de la cola Llegada de un cliente Cola Servicio
Fuente de entrada
Mecanismo de servicio Sistema de la cola
Figura 2- Elementos existentes en un modelo de colas
19
Capítulo 1: Introducción
Jorge L. Vega Valle
1.2.5- Terminología y notación 1.2.5.1- Terminología
A menos que se establezca otra cosa, y tomando como referencia [Cao-02], se utilizará la siguiente terminología estándar: N(t): Denota el número de clientes en el sistema en el instante t. N(t) es un
proceso estocástico en tiempo continuo y con espacios de estados discreto. Nq(t): Representa el número de clientes en la cola en el instante t. Pn(t): Es la probabilidad de que, en el instante t, se encuentren n clientes en
el sistema. A estos efectos se supone conocido el número de clientes en el instante cero (usualmente dicho número es cero). s: Denota el número de servidores del mecanismo de servicio. λn:
Representa el número medio de llegadas de clientes al sistema, por
unidad de tiempo, cuando ya hayn clientes en él. También se denomina tasa de llegadas (que se correspondería con la tasa de nacimientos si N(t) es un proceso de nacimiento y muerte). Cuando las tasas de llegada no dependen den (es decir todos los λn son constantes) suele denotarse por λ dicho valor constante. µn: Es el número medio de clientes a los que se les completa el servicio,
por unidad de tiempo, cuando hay n clientes en el sistema. Es frecuente referirse a los µn como tasas de complección de servicio (o, simplemente, tasas de servicio). Si todos los servidores tienen la misma distribución del tiempo de servicio, suele denotarse por µ el número medio de clientes que puede atender
20
Capítulo 1: Introducción
Jorge L. Vega Valle
cada servidor por unidad de tiempo. Como consecuencia se tiene queµn= n • µ si n = 1, 2,..., s y µn = s • µ para n ≥ s. ρ:
Es la llamada constante de utilización del sistema o intensidad de
tráfico. Se define, como
ρ=
λ sµ
Cuando los λn son constantes y todos los servidores tienen la misma distribución de tiempo de servicio, λ es el número medio de clientes que entran en el sistema y sµ es el número medio de clientes a los que pueden dar servicio los s servidores cuando todos están ocupados. En estas condiciones, ρ representa la fracción de recursos del sistema que es consumida por los clientes. Así, intuitivamente, parece necesario que se cumpla, en estos casos, que ρ < 1 y además cuanto más cercano a 1 que sea su valor, más tráfico ha de soportar el sistema (o menos tiempo libre tendrán los servidores, o más espera habrán de sufrir los clientes, como se quiera expresar). Aunque es evidente que ρ no tiene unidades, es habitual medir la intensidad de tráfico en Erlangs, en honor a los trabajos pioneros de Erlang en la teoría de colas.
1.2.5.2- Conceptos básicos Los siguientes conceptos, como se puede ver en [Cao-02], son de utilidad para analizar las características y el comportamiento de un modelo de colas estacionario:
N: Es la variable aleatoria que contabiliza el número de clientes en el sistema. Nq: Denota la variable aleatoria número de clientes en la cola.
21
Capítulo 1: Introducción
Jorge L. Vega Valle
pn: Es la probabilidad de que se encuentren n clientes en el sistema (n = 0,
1,…). L: Representa el número medio de clientes en el sistema, es decir L = E(N). Lq: Que no es más que el número medio de clientes en la cola, o lo que es lo
mismo, Lq = E(Nq). : Es la variable aleatoria que describe el tiempo que un cliente pasa en el sistema o también llamado tiempo de espera en el sistema (incluyendo el tiempo de servicio) para cada cliente.
q
: Representa el tiempo que un cliente espera en la cola.
W: Es el tiempo medio que un cliente está en el sistema. En términos matemáticos, W = E( ).
Wq: Denota el tiempo medio de espera en la cola para un cliente genérico.
Matemáticamente, Wq = E( q).
1.2.5.3- Notación de Kendall Para clasificar los posible tipos de sistemas de colas debemos especificar las características que determinan los elementos que lo componen. Así, Kendall (ver [Bos-02], [Cao-02] o [Hil-97]) introdujo en 1953 la notación A/B/s para indicar que la distribución del tiempo entre llegadas es de del tipo A, que B es la distribución del tiempo de servicio y que s es el número de servidores. Posteriormente esta notación
22
Capítulo 1: Introducción
Jorge L. Vega Valle
se extendió dando lugar a la más habitual en nuestros días, consistente en designar el sistema de una cola con la nomenclatura A/B/s/K/H/Z, donde: A es la distribución del tiempo entre llegadas. Algunas de las abreviaturas
más usadas para las distribuciones entre llegadas son: M (exponencial), D (determinística), Ek (Erlang con segundo parámetrok), U (uniforme), Γ (gamma) o G (distribución genérica), entre otras. B es la distribución del tiempo de servicio. Se usan las mismas abreviaturas
que las mencionadas paraA. s es el número de servidores del sistema. Puede ser un número entero positivo
(s = 1, 2,…) o bien s = ∞ . K es la capacidad de la cola (o longitud máxima de la misma). También K
puede ser un número entero mayor o igual que cero, o bien K = ∞ , si no hay límite para la cola. El valor de K puede omitirse, tomándose por defecto K = ∞ . H es el tamaño de la población potencial. También puede ser finito o infinito.
Este último valor es el que se toma por defecto cuando se omite su valor. Z es la disciplina en la cola. Algunas abreviaturas para Z son FIFO, LIFO,
RSS, PR (disciplina con prioridades) o GD (disciplina general). Su valor por defecto (en caso de omitirse Z) es FIFO. Así, por ejemplo, la notación M/D/2/∞ / ∞ /FIFO indica que se trata del sistema de una cola con tiempo entre llegadas exponenciales, tiempo de servicio determinístico (i.e. siempre se tarda el mismo tiempo en darle servicio a cada cliente), hay 2 servidores en el mecanismo de servicio, no existe límite para el número de clientes que pueden estar en la cola de espera, la población potencial se supone con infinitos
23
Capítulo 1: Introducción
Jorge L. Vega Valle
clientes y los clientes son atendidos según una disciplina FIFO. Como los tres últimos valores (∞ ,
∞
y FIFO) son precisamente los asignados por defecto, la
notación anterior podría abreviarse como M/D/2.
1.2.6- Redes de colas Se presentan aquí los modelos básicos de redes de colas abiertas y cerradas con distribución del tiempo de servicio y distribución del tiempo entre llegadas (si es el caso) exponencial. Además se impondrá la restricción de que los clientes que salen servidos de una de las colas que compone la red se mueven instantáneamente y con ciertas probabilidades prefijadas a cualquier otra posible cola de la red. Estos modelos dan lugar a las llamadas redes de Jackson (abiertas y cerradas).
1.2.6.1- Introducción a las redes de colas Una red de colas no es más que una red en la que cada nodo está constituido por el sistema de una cola. Se trata, por tanto, de un grafo orientado en el que se pueden producir transiciones de clientes que salen servidos de un nodo (que es una cola) hacia otro nodo. La forma más habitual (aunque no la única) para modelizar el modo en que los clientes servidos en un nodo se dirigen a otro es considerando que lo hacen de acuerdo a una distribución de probabilidad discreta. Al igual que en los modelos de una única cola, en las redes de colas abiertas también pueden producirse llegadas de clientes desde fuera del sistema (desde fuera de la red, en este caso) y salidas de clientes servidos hacia fuera de la red. A diferencia de aquél caso, en las redes de colas sí tiene sentido el plantear situaciones en las que no hay llegadas de clientes desde fuera de la red ni salidas hacia fuera de
24
Capítulo 1: Introducción
Jorge L. Vega Valle
la red. Esto da lugar a las llamadas redes cerradas y tienen la peculiaridad que el número total de clientes en la red es fijo, y lo único que desconocemos es dónde se encuentran (en qué nodos concretos) y en qué estado de servicio se hallan. Por su parte, las redes abiertas son aquellas en que sí se producen llegadas de clientes y salidas hacia fuera de la red. Como se comentaba con anterioridad, denotando por 1, 2,…, K, los nodos que forman la red (es decir las etiquetas con las que denotamos a cada cola) la manera más frecuente de modelizar las transiciones de clientes consiste en suponer que cuando un cliente sale servido de la cola del nodo i (siendo i ε{1, 2,…, K}) se desplaza instantáneamente al sistema de la cola de cualquier otro nodo j
ε
{1,
2,…, K}, con probabilidad pij. Evidentemente, en las colas abiertas también es posible que desde algunos nodos se pueda abandonar la red. Denotando con el índice 0 el exterior de la red, la probabilidad de que un cliente abandone la red cuando sale servido del nodo i se denotará por pi0 y puede calcularse a partir de las anteriores mediante pi 0 = 1 −
K
∑1 p
ij
j=
Este tipo de esquema de transición de clientes se denotará esquema de transiciones instantáneas aleatorias de clientes.
En una situación como la anterior, las probabilidades de transición de clientes de unos nodos a otros se pueden expresar de forma matricial:
25
Capítulo 1: Introducción
p11 p 21 M P= p i1 M K p 1
Jorge L. Vega Valle
p12
L
p1 j
L
p1K
p 22
L
p2 j
L
p2 K
M
O
M
O
pi 2
L
p ij
L
M
O
M
O
pK 2
L
p Kj
L
p iK M p KK M
Esta es la llamada matriz de transición de la red. Normalmente denotaremos por λi la tasa de entrada desde fuera del sistema a la cola del nodo i ε {1, 2,…, K}. Asimismo µi denotará la tasa de servicio de cada uno de los servidores del subsistema del nodo i. De esta forma, bajo la hipótesis de que los tiempos entre llegadas desde fuera del sistema sean de distribución exponencial de parámetro λT (que denota la tasa de entrada al sistema, λT = Σ λi ) y denotando por p0i la probabilidad de que cuando un cliente entra al sistema lo haga a través del nodo i, se tiene que la distribución del tiempo entre dos entradas consecutivas de clientes
al subsistema del nodo i es también exponencial y con parámetro λi = λT • p0i. Obviamente los valores de λi también pueden interpretarse como el número medio de clientes que entran al sistema, por el nodo i, por unidad de tiempo. El número efectivo de llegadas de clientes al subsistema del nodo i (sean procedentes de fuera
del sistema o de otro nodo del mismo) se denotará con la letra lambda mayúscula Λi.
1.2.6.2- Redes de Jackson abiertas Una red de Jackson abierta no es más que una red de colas abierta (es decir, en la es posible la llegada de clientes desde fuera de la red y la salida de clientes a fuera de la red) que verifica las tres propiedades siguientes:
26
Capítulo 1: Introducción
Jorge L. Vega Valle
Cada nodo i = 1, 2,…, K tiene un mecanismo de servicio consistente en si servidores con tiempo de servicio de idéntica distribución exponencial de parámetro µi, de tal forma que las distribuciones de los tiempos de servicio de cada nodo son iguales en todos los servidores de un nodo (pudiendo ser distintos entre los diferentes nodos de la red). Los clientes que llegan al nodo i desde fuera del sistema lo hacen según un proceso de Poisson de intensidad λi. Esto equivale a decir que los tiempos entre
dos llegadas de clientes consecutivos desde fuera del sistema al nodo i, siguen una distribución exponencial de dicho parámetro. El flujo de clientes sigue el esquema de transiciones instantáneas aleatorias. Es decir, cada cliente que sale servido del mecanismo de servicio del nodo i va instantáneamente a cualquier otro nodo j con probabilidad pij, o bien sale del K
sistema (lo cual ocurre con la siguiente probabilidad: pi 0 = 1 − ∑ p).ij j =1
En la siguiente figura podemos ver un ejemplo de red de colas abierta:
27
Capítulo 1: Introducción
Jorge L. Vega Valle
p11
p22
p12
λ1 NODO 1
p20
NODO 2
p21
λ2
p10 p32 p13
p23 p31
EXTERIOR
EXTERIOR
λ3 NODO 3 p30 p33
Figura 3- Red de Jackson abierta con K = 3
1.2.6.3- Redes de Jackson cerradas Una red de Jackson cerrada es una red de colas cerrada conK nodos o subsistemas, en la cual cada nodo i = 1, 2,…, K tiene si servidores en su mecanismo de servicio, siendo todos los del nodo i con tiempo de servicio de distribución exponencial de parámetro µi. A diferencia de las redes abiertas, no es posible ni la entrada ni salida de clientes hacia el exterior, con lo que, resulta indispensable especificar en número de clientes dentro de la red, N, que permanecerá constante siempre. Por este motivo, LT = N y cantidades como WT o Wq,T, carecen de sentido. Lo realmente importante
aquí es determinar las probabilidades de que haya ni clientes en el nodo i para i = 1, 2,…, K, que se denotarán por p n1 ,n2 ,...n . Obviamente éstas probabilidades son distintas de cero sólo si n1 + n2 + ••• + nK = N. K
28
Capítulo 1: Introducción
Jorge L. Vega Valle
Como no hay entradas de clientes, las probabilidades de transiciones entre nodos deben verificar lo siguiente: k
∑1 p
ij
=1
(siendo i un nodo cualquiera de la red de colas)
j=
Tal y como se comenta en [Bos-02], una red de colas cerrada puede parecer inusual a primera vista debido a que no permite entrada o salida de clientes de la red. Un sistema típico de este tipo podría ser uno en el cual haya una gran cantidad (infinita) de clientes intentando entrar al sistema de forma continua, el número de clientes a los que se les permite la entrada tiene un determinado valor fijo y un cliente entra en el sistema de forma inmediata, siempre que se complete la secuencia de servicios de otro cliente (por ejemplo, un sistema de computación de tiempo compartido). En la siguiente figura podemos ver un ejemplo de red de colas cerrada:
29
Capítulo 1: Introducción
Jorge L. Vega Valle
p11
p22
p12
NODO 1
NODO 2
p21 p32
p13
p23 p31
NODO 3
p33
Figura 4- Red de Jackson cerrada con K = 3
30
Capítulo 1: Introducción
Jorge L. Vega Valle
1.3- Introducción a la Simulación Desde hace mucho tiempo la simulación ha sido una herramienta importante en una gran cantidad de ámbitos. Por ejemplo, la simulación del vuelo de un avión en un túnel de viento es una práctica normal cuando se diseña un nuevo avión. En teoría se podrían usar las leyes de la física para obtener la misma información sobre los cambios que sufre el avión si se modifican los parámetros, pero en el sentido práctico, el análisis sería muy complicado. Una alternativa sería construir aviones reales para cada uno de los diseños y probarlos en vuelos reales, pero esto sería demasiado costoso (a la vez que peligroso). Por lo tanto, después de realizar un análisis teórico preliminar para desarrollar un diseño básico, la herramienta más viable para experimentar con los diseños específicos es la simulación del vuelo en un túnel de viento. Esta simulación significa imitar el desempeño de un avión real en un medio controlado con el fin de estimar cuál sería el desempeño real. Después de desarrollar un diseño detallado se puede construir un modelo prototipo y probarse en un vuelo real para dar los últimos detalles del diseño final.
1.3.1- Conceptos básicos La simulación es la técnica que consiste en realizar experimentos de muestreo sobre el modelo de un sistema. Un modelo no es más que un conjunto de variables junto con ecuaciones matemáticas que las relacionan y restricciones sobre dichas variables. La modelización es una etapa presente en la mayor parte de los trabajos de investigación (especialmente en las ciencias experimentales). En muchas ocasiones, la realidad es bastante compleja como para ser estudiada directamente y es preferible la formulación de un modelo que contenga las variables más relevantes que aparecen en el fenómeno en estudio y las relaciones más importantes entre ellas.
31
Capítulo 1: Introducción
Jorge L. Vega Valle
Frecuentemente, la resolución de los problemas que se pretenden abordar puede realizarse por procedimientos analíticos sobre el modelo construido (normalmente mediante el uso de herramientas matemáticas como las de resolución de ecuaciones ordinarias o de ecuaciones diferenciales, el cálculo de probabilidades, etc.). En otras circunstancias dicha resolución analítica no es posible (o es tremendamente complicada o costosa) y es preferible una aproximación de la solución mediante simulación.
1.3.2- Experimentación real y simulación Como se comenta en [Cao-02], la experimentación directa sobre la realidad puede tener muchos inconvenientes: un coste muy alto gran lentitud en ocasiones las pruebas son destructivas a veces no es ética (experimentación sobre seres humanos) puede resultar imposible (un acontecimiento futuro) Razones como esas (y algunas otras) pueden indicar la ventaja de trabajar con un modelo del sistema real. La estadística es precisamente la ciencia que se preocupa de cómo estimar los parámetros y contrastar la validez de un modelo a partir de los datos observados del sistema real que se pretende modelizar.
32
Capítulo 1: Introducción
Jorge L. Vega Valle
1.3.3- Tipos de simulación 1.3.3.1- Simulación estática y dinámica
La simulación se dice estática si en el modelo no juega ningún papel el transcurso del tiempo mientras que es dinámica si el tiempo es una de las variables importantes del modelo. En la simulación estática resulta muy sencillo comparar distintas estrategias ante las mismas condiciones del azar, mientras que esto es más complicado en la simulación dinámica, exigiendo un trabajo mayor de planificación. Además, el coste computacional de la simulación estática es bastante más moderado. La simulación estática se usa muy frecuentemente por los estadísticos para comprobar el comportamiento comparativo de diversos métodos estadísticos alternativos para tamaños muestrales finitos (complementando los estudios teóricos, casi siempre asintóticos). En la simulación dinámica, normalmente se trata de ir analizando los distintos estados por los que va pasando un sistema que evoluciona en el tiempo. Esto provoca, en general, un mayor coste computacional y problemas de estabilización y dependencia. Existen dos grandes tipos de simulación dinámica: la simulación continua, en la que se supone que el sistema cambia de estado constantemente y la simulación discreta, para la cual los cambios se producen en ciertos instantes de tiempo singulares. La razón de sus nombres viene de que en el primer caso el conjunto de estados es continuo, mientras que en el segundo es discreto. Dentro de la simulación discreta distinguiremos la simulación por eventos y la simulación por cuantos.
33
Capítulo 1: Introducción
Jorge L. Vega Valle
1.3.3.2- Simulación por eventos y por cuantos Con el nombre de simulación por eventos, o asíncrona, designamos el tipo de simulación dinámica discreta en la cual se controla la variable tiempo moviéndola hasta la ocurrencia del siguiente suceso (o evento). Esto implica la necesidad de controlar minuciosamente cuál es dicho próximo suceso: saber cuáles son los posibles sucesos en un futuro inmediato y cuál de ellos es el más inmediato. Resumen de la simulación por eventos: [Hil-97]
Se avanza el tiempo hasta el momento en que ocurre el siguiente evento de cualquier tipo. Se actualiza el sistema determinando su nuevo estado, que es el resultado de este evento y generado aleatoriamente (si no se generó antes) el tiempo hasta la siguiente ocurrencia de un evento de cualquier tipo que pueda ocurrir estando en este estado. También se registra la información deseada sobre el comportamiento del sistema. (volver al paso anterior). Ejemplo de simulación por eventos: [Bos-02]
Supongamos una red de colas abierta con dos nodos. Los eventos que debemos considerar son los siguientes: Llegadas del exterior al nodo 1 Llegadas del exterior al nodo 2 Llegadas al nodo 1 procedentes del nodo 1 Llegadas al nodo 1 procedentes del nodo 2 Llegadas al nodo 2 procedentes del nodo 1
34
Capítulo 1: Introducción
Jorge L. Vega Valle
Llegadas al nodo 2 procedentes del nodo 2 Finalización del servicio (salida) en el nodo 1 Finalización del servicio (salida) en el nodo 2
Como se comenta en [Cao-02], la simulación por cuantos, o asíncrona, responde a una filosofía totalmente diferente. Se trata de examinar el sistema (que evoluciona en el tiempo) dejando pasar pequeños intervalos de tiempo de longitud δ, fija, (llamada cuanto) en los cuales se supone que, a lo sumo, un sólo suceso puede producirse. Resumen de la simulación por cuantos: [Hil-97]
Se avanza el tiempo una cantidad fija pequeña.
Se actualiza el sistema determinando los eventos que ocurrieron durante ese lapso y el estado del sistema que resulta. También se registra la información deseada sobre el comportamiento del sistema (volver al paso anterior). En general, la simulación por eventos es exacta y de más difícil implementación, pero de mucha más rápida ejecución que la simulación por cuantos. Sin embargo esta última es muchas veces la única posibilidad factible en la simulación dinámica continua.
35
Capítulo 1: Introducción
Jorge L. Vega Valle
En el siguiente esquema podemos ver un resumen de los tipos de simulación:
Estática
Simulación
{
Dinámica
{
Continua
Discreta
{
Por Eventos
Por Cuantos
Figura 5- Resumen de tipos de simulación
1.3.4- Problemas de estabilización y dependencia Ambas cuestiones suelen plantearse en la simulación dinámica. Los problemas de estabilización están relacionados con el hecho de que, en ocasiones, el sistema evoluciona en el tiempo de tal forma que tiene una distribución estacionaria que se supone de partida pero que puede ser muy sensible a las condiciones iniciales con las que se comience la simulación. En tal caso resulta conveniente el transcurso de cierto período de tiempo (denominado período de estabilización) durante el cual los resultados obtenidos para las variables de salida son ignorados y cuyo único objeto es conseguir que se estabilice la distribución de probabilidad. Los problemas de dependencia son aquellos derivados del hecho de que frecuentemente (de nuevo en modelos de simulación dinámica) las distintas variables
36
Capítulo 1: Introducción
Jorge L. Vega Valle
de salida de la simulación son dependientes. Esto afecta fundamentalmente a la precisión de los estimadores construidos con observaciones de las mismas. Una forma de atenuar este efecto sería considerar observaciones de las mismas en instantes temporalmente lejanos (donde se supone que la dependencia es mucho más débil). En ocasiones, más que atenuar este efecto se trata de estimar la precisión del estimador resultante. Obviamente, para ello ha de tenerse en cuenta la dependencia.
1.3.5- Comportamiento de la simulación En la siguiente figura, se puede observar el comportamiento típico de la simulación (descrito con más detalle en [Bos-02]) de un parámetro Z durante una simulación que fue llevada a cabo en un intervalo de tiempo (0, T). El parámetro Z se encuentra inicialmente en un estado transitorio, no necesariamente cercano a su verdadero valor, pero en un momento dado alcanza valores ya muy cercanos a la cantidad de interés. En la figura se pueden observar dos intervalos: el primero de ellos, llamado intervalo de transición transcurre cuando t se encuentra entre 0 y TZ, mientras que el segundo de ellos, en el cual la simulación ha alcanzado las condiciones de equilibrio, transcurre cuando t se encuentra entre TZ y T.
37
Capítulo 1: Introducción
Jorge L. Vega Valle
Estado de Equilibrio Estado Transitorio Simulación Variable, Z
0
T
TZ
Tiempo
Figura 6- Comportamiento de la simulación
1.3.6- Aplicaciones comunes y uso de la simulación Se han realizado numerosas aplicaciones de simulación en una gran variedad de contextos. Estos son sólo algunos ejemplos (podemos encontrar más en [Hil-97]) que ilustran la gran versatilidad de esta técnica: Simulación de las operaciones en un aeropuerto por una compañía aérea para probar cambios en las políticas y prácticas de la empresa. Simulación de la economía para predecir el efecto de las decisiones de las políticas económicas.
38
Capítulo 1: Introducción
Jorge L. Vega Valle
Simulación de batallas militares de gran escala para evaluar los sistemas de armamento ofensivo y defensivo. Simulación de un sistema de comunicaciones telefónicas para determinar la capacidad requerida de las respectivas componentes para proporcionar un servicio satisfactorio al nivel más económico. Simulación de la operación de la cuenca hidráulica de un río para determinar la mejor configuración de presas, plantas de energía y sistemas de irrigación que proporcione el nivel deseado de control de avenidas y desarrollo del recurso. La simulación es un método de análisis científico ampliamente usado. Un estudio comentado por [Wat-89], de The Institute of Managemet Sciences (TIMS) y de la Operactions Research Society of America (ORSA) averiguó que la simulación es una de las técnicas científicas que más se usa, solo detrás de las técnicas de análisis estadístico (inferencia, teoría de decisiones,...) y las técnicas de análisis económico. Otro estudio de los miembros de TIMS / ORSA comprobó si las organizaciones usaban técnicas de simulación, y, si efectivamente era así, que resultados habían obtenido. Los resultados fueron que el 89% de las organizaciones utilizaban técnicas de simulación y que el 85% de las mismas estaban contentas con los resultados obtenidos.
39
Capítulo 2: Aspectos de implementación “No queda sino batirnos” (El personaje que protagoniza Francisco de Quevedo en la novela “El Capitán Alatriste” de Arturo Pérez-Reverte)
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.1- Características generales de la aplicación Como ya se ha comentado anteriormente, la aplicación desarrollada es capaz de resolver problemas de colas de forma analítica y mediante simulación. Por ello, hemos decidido llamar a la aplicación con el siguiente nombre: AQUAS, cuyas siglas significan Application for solving QUeuing problems Analytically and using Simulation.
2.1.1- Filosofía de la implementación Como ya se ha comentado en anteriores apartados (ver 1.1.3), este proyecto ha sido desarrollado
utilizando
MATLAB®
como
herramienta
fundamental.
La
programación de la aplicación ha sido realizada utilizando tanto ficheros de comandos (scripts), como funciones, mientras que el control de los distintos eventos se ha desarrollado empleando “callbacks1 ”. Es importante reseñar que, a pesar de que la palabra “diseño” está incluida en el titulo del proyecto, no tiene la connotación típica en informática (ya que esta suele ir asociada a la programación orientada a objetos), sino que más bien se refiere al desarrollo y elección de las características estéticas de la aplicación. Por otra parte, para la realización de la interfaz se ha usado la GUI2 de MATLAB y las características que proporciona. También se han empleado algunas funciones específicas de la librería estadística (ver 2.1.4 y [2]).
1 2
Para una mayor aclaración sobre estos conceptos, véase [Red-98] o [8] GUI: Graphical User Interface. Ver [Mar-99] o [3]
41
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.1.2- Ciclo de vida empleado Hay que diferenciar entre el ciclo de vida (definición de [Pre-97]) usado para la parte analítica y para la parte de simulación. En la parte analítica, el análisis de cada modelo es muy sencillo, ya que los algoritmos de resolución están previamente establecidos, cosa que no ocurre en la parte de simulación, lo que implica que el análisis sea más complejo, ya que hay que tener en cuenta tanto los posibles eventos (ver 1.3.3.2) que pueden ocurrir como la influencia de estos eventos en el cálculo de los distintos parámetros de salida. En las siguientes figuras podemos observar con detalle el ciclo de vida usado para la parte analítica y para la parte de simulación.
DETERMINAR LOS PARAMETROS DE ENTRADA ERROR
DETERMINAR LOS PARAMETROS DE SALIDA IMPLEMENTAR EL MODELO DETERMINAR COMO CALCULAR LOS PARAMETROS DE SALIDA
ANALIZAR LOS POSIBLES RIESGOS
ANALISIS DE UN DETERMINADO MODELO
Figura 7- Ciclo de vida para la parte de resolución analítica
42
PRUEBAS SW DEL MODELO IMPLEMENTADO
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
DETERMINAR LOS PARAMETROS DE ENTRADA ERROR
DETERMINAR LOS PARAMETROS DE SALIDA
ANALIZAR LA IMPLEMENTACION DEL ALGORITMO
ANALIZAR LOS POSIBLES RIESGOS
IMPLEMENTAR EL MODELO
ESTUDIAR LOS POSIBLES EVENTOS QUE PUEDEN OCURRIR (llegada de cliente, salida, transición)
PRUEBAS SW DEL MODELO IMPLEMENTADO (CONTRASTE CON RESOLUCION ANALITICA)
ANALIZAR LA INFLUENCIA DEL EVENTO EN LA VARIACION DE LOS PARAMETROS
ANALISIS DE UN DETERMINADO MODELO
Figura 8- Ciclo de vida para la parte de simulación
2.1.3- Parámetros de entrada y salida La aplicación da respuesta a los distintos modelos de colas. Cuando la hipótesis de exponencialidad es cierta, se proporciona una respuesta analítica exacta al modelo. Por otra parte, si no podemos suponer dicha hipótesis, deberemos recurrir a la simulación para poder tener una respuesta con las características del modelo. Tanto para la resolución analítica como en la resolución por simulación los parámetros de entrada son los siguientes: Número de servidores (s):
Sólo para el caso de los modelos M/M/s, M/M/s/K, M/M/s/∞ /H y M/M/s/∞ /H con Y repuestos en resolución analítica y los modelos G/G/s,
43
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
G/G/s/K, G/G/s/∞ /H y G/G/s/∞ /H con Y repuestos en el caso de resolución
por simulación. En el caso de que estemos en una red de colas debemos introducir el número de servidores para cada uno de los nodos de la red. Tamaño de la cola (K): Sólo para el caso de los modelos M/M/1/K y M/M/s/K en resolución analítica y los modelos G/G/1/K y G/G/s/K en el caso de
resolución por simulación. Población potencial (H): Sólo para el caso de los modelos M/M/1/∞ /H, M/M/s/∞ /H y M/M/s/∞ /H con Y repuestos en resolución analítica y los modelos G/G/1/∞ /H, G/G/s/∞ /H y G/G/s/∞ /H con Y repuestos en el caso de resolución
por simulación.
Número de repuestos (Y): Sólo para el caso del modelo M/M/s/∞ /H con Y repuestos en resolución analítica y el modelo G/G/s/∞ /H con Y repuestos en el
caso de resolución por simulación. Probabilidades de transición entre nodos (pij): En el caso de las redes de
colas abiertas o cerradas. Número de nodos que forman la red (K): En el caso de que nos
encontremos en una red de colas Número de clientes de la red cerrada (N)
En el caso de que el usuario seleccione la opción deresolución analítica tendrá además los siguientes parámetros de entrada:
44
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Tasa de llegada (λ):
En las redes de colas abiertas debemos introducir la tasa de llegada para cada uno de los nodos. En las redes de colas cerradas no tiene sentido el incluir este parámetro. Tasa de servicio (µ): En las redes de colas debemos de incluir la tasa de
servicio para cada uno de los nodos. Mientras que si elige la resolución por simulación tiene, a mayores, los siguientes parámetros de entrada: Tipo de distribución a la entrada y la salida: En el caso de las redes de
colas no tiene sentido la distribución de entrada. Parámetros asociados a las distintas distribuciones elegidas. Parámetro de estabilización. Número de clientes de la simulación ó Numero de transiciones de clientes.
Por otra parte, los parámetros de salida son los siguientes: L: Número medio de clientes en el sistema. En el caso de una red de Jackson
abierta se muestra el número medio de clientes en cada nodo y en toda la red. En el caso de una red de Jackson cerrada se muestra el número medio de clientes en cada nodo de la red.
45
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Lq: Número medio de clientes en cola. En el caso de una red de Jackson
abierta se muestra el número medio de clientes en la cola en cada nodo y en toda la red. En el caso de una red de Jackson cerrada se muestra el número medio de clientes en la cola en cada nodo de la red. W: Tiempo medio que un cliente está en el sistema. En el caso de una red de
Jackson abierta se muestra el tiempo medio que un cliente está en el sistema en cada nodo y en toda la red. En el caso de una red de Jackson cerrada se muestra el tiempo medio que un cliente está en el sistema en cada nodo de la red. Wq: Tiempo medio de espera en la cola para un cliente. En el caso de una red
de Jackson abierta se muestra el tiempo medio de espera en la cola para un cliente en cada nodo y en toda la red. En el caso de una red de Jackson cerrada se muestra el tiempo medio de espera en la cola para un cliente en cada nodo de la red. Intensidad del modelo: Refleja la intensidad de tráfico (también conocida
como ρ ) y toma valores entre 0 y 1. Es necesario que ρ sea inferior a 1, para que el modelo sea estacionario. No se muestra en las redes de colas. Eficiencia del modelo: Equivale al cociente entre el tiempo medio que un
cliente está en el sistema (W) y el tiempo medio de servicio (Ws), por lo que su valor siempre es mayor o igual que 1. No se muestra en las redes de colas. Valor de W(t) en un instante determinado: No se muestra ni en las redes de
colas ni en el caso de resolución por simulación. Valor de Wq(t) en un instante determinado: No se muestra ni en las redes
de colas ni en el caso de resolución por simulación.
46
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Gráfica con los valores de W(t) y Wq(t) en unos instantes determinados:
No se muestra ni en las redes de colas ni en el caso de resolución por simulación. p(n): Probabilidad de que hayan clientes en el sistema.
En el caso de las redes abiertas y resolución analítica podemos hallar la probabilidad de que haya un número determinado de clientes en cada uno de los nodos de la red. En el caso de las redes cerradas y resolución analítica tenemos la posibilidad de hallar la probabilidad de que haya n clientes en un nodo determinado de la red, o hallar la probabilidad de que haya un número determinado de clientes en cada uno de los nodos que forman la red. En el caso de las redes abiertas o cerradas y resolución por simulación tenemos la posibilidad de hallar la probabilidad de que haya un número determinado de clientes en un nodo de la red. q(n): Probabilidad de que haya n clientes en el sistema justo cuando se
produce una llega. No se muestra ni en los modelos M/M/1, M/M/s en el caso de resolución analítica, ni en los modelos G/G/1, G/G/s en el caso de resolución por simulación, ni en las redes de colas. Gráfica con los valores de p(n) cuando n toma unos valores determinados. No se muestra en el caso de las redes abiertas, ni en la simulación de redes cerradas.
47
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.1.4- Distribuciones de probabilidad elegidas Para resolver la parte de simulación hemos elegido un conjunto de distribuciones de probabilidad3 representativas. El usuario puede variar tanto la distribución de llegada de clientes como la distribución de servicio. Las distribuciones elegidas son las siguientes (se explicarán con más detalle en el Apéndice A): Exponencial Uniforme Determinista Gamma Beta Lognormal Normal De Weibull Para generar números aleatorios de dichas distribuciones se han usado las funciones (ver [2]) que implementa MATLAB. En la siguiente tabla podemos observar dichas funciones. Nótese que la definición de MATLAB para algunas de las funciones, como por ejemplo la gamma o la de Weibull, es ligeramente distinta a la nuestra:
3
Podemos observar su definición en [Cao-98] o [Dev-86]
48
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
DISTRIBUCION
PARAMETROS
FUNCION MATLAB
EXPONENCIAL
λoµ
exprnd(1/ λ) o exprnd(1/ µ)
UNIFORME
a, b
unifrnd(a, b)
DETERMINISTA
d
No hay. Es trivial.
GAMMA
a, p
gamrnd(p, 1/a)
BETA
p, q, k
k • betarnd(a, b)
LOGNORMAL
µ, σ
lognrnd(µ, σ)
NORMAL
µ, σ
normrnd(µ, σ)
DE WEIBULL
λ, α
weibrnd(λα, α)
Es importante resaltar que en ningún momento se establece un valor para la semilla. Cada vez que se llame a la función de generación de números de una determinada distribución, se utilizará el reloj del sistema como semilla, lo que implica que si se hacen varias llamadas consecutivas a la función se obtendrán valores no necesariamente iguales.
49
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.2- Resolución analítica 2.2.1- Características generales
Para resolver de forma analítica los distintos modelos de colas se han implementado las soluciones planteadas por los distintos libros (véase [All-90], [Bos-02], [Cao-02], [Gro-98] o [Hil-97]). La implementación de dichas fórmulas se realizó de forma eficiente, utilizando, en la medida de lo posible, la versión de las fórmulas cuyo coste computacional era más bajo y evitando la implementación directa de las fórmulas, de escasa eficiencia computacional. En el siguiente apartado se explica de forma detallada las fórmulas que resuelven cada uno de los modelos de colas (modelos de la forma M/M/- y redes de Jackson).
2.2.2- Resolución de los modelos 2.2.2.1- Modelo M/M/14 Intensidad = ρ =
L=
4
λ
Eficiencia =
µ
λ
W W = =W ⋅µ 1 W − Wq µ
Lq =
µ −λ
λ2 µ ⋅ (µ − λ )
La fórmula de la eficiencia es siempre igual, por lo que, se obviará en los siguiente modelos
50
Capítulo 2:Aspectos de implementación
W=
Jorge L. Vega Valle
1
Wq =
µ −λ
W (t ) = 1 − e −(µ− λ) ⋅ t
si t ≥ 0
Wq (t ) = 1 −
λ
λ µ ⋅ (µ − λ )
⋅ e −(µ− λ) ⋅ t si t ≥ 0
µ
p n = (1 − ρ ) ⋅ ρ n si n = 0, 1, ...
2.2.2.2- Modelo M/M/s Intensidad = ρ =
c0 = 1
λ s⋅µ
cn = cn−1 ⋅ c≥ s = c s−1 ⋅
p0 =
λ
⋅
si n = 1, 2,..., s - 1
n µ λ s⋅µ −λ
1 s −1
∑0 c
n
+ c≥ s
n=
p n = c n ⋅ p0 pn = ρ ⋅ p n−1
si n = 1, 2,..., s - 1 si n ≥ s
Lq = c ≥ s ⋅ λ ⋅ p 0 s⋅µ −λ
L = λ ⋅W
51
Capítulo 2:Aspectos de implementación
W = Wq +
Jorge L. Vega Valle
1
Wq =
µ
Lq λ
Wq (t ) = 1 − c≥ s ⋅ p0 ⋅ e − ( s⋅µ −λ )⋅t si t ≥ 0
W (t ) = 1 − (1 + c≥ s ⋅ p0 ⋅ t ⋅ µ ) ⋅ e − µ⋅t W (t ) = 1 +
λ − s ⋅ µ + µ ⋅ Wq (0)
s⋅µ −λ − µ
si t ≥ 0 y
⋅ e − µ⋅t +
λ µ
= s −1
c≥ s ⋅ µ ⋅ p0 ⋅ e −( s⋅µ −λ )⋅t s⋅µ −λ −µ
si t ≥ 0 y
λ µ
≠ s −1
2.2.2.3- Modelo M/M/1/K
Intensidad = ρ =
ρ=
λ
λ=
µ
λ
λ=
µ
pn = pn =
1 K+2 ρ −1
qn =
−1 pn
1 − pK +1
( K + 2)
λ ⋅ ( ρ K +1 − 1) ρ K +2 − 1
si n = 0, 1,..., K + 1 y ρ = 1
K +2
ρ
λ ⋅ ( K + 1)
⋅ ρn
si n = 0, 1,..., K + 1 y ρ ≠ 1
si n = 0, 1,..., K
52
si ρ = 1
si ρ ≠ 1
Capítulo 2:Aspectos de implementación
K +1
L=
L=
1− ρ
W =
Lq = λ ⋅ W q
si ρ = 1
2 ρ
Jorge L. Vega Valle
−
( K + 2) ⋅ ρ K + 2 1 − ρ K +2
si ρ ≠ 1
L
Wq = W −
λ
Cálculo de W(t):
Cálculo de Wq(t):
A = 1, S = 1, B = q0
A = 1, S = 1, B = q1
Desde n = 1 hasta K repetir
Desde n = 2 hasta K repetir
A = A • (µ • t) / n
A = A • (µ • t) / n
S=S+A
S=S+A
B = B + qn • S Devolver 1 – B • e –µ • t
B = B + qn • S Devolver 1 – B • e –µ • t
53
1 µ
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.2.2.4- Modelo M/M/s/K5 c0 = 1
Intensidad = ρ = ρ ⋅ (1 − p K + s )
cn = cn−1 ⋅ ρ=
λ n⋅µ
si n = 1, 2,..., s - 1
c≥ s = c s −1 ⋅ ( K + 1)
λ
s⋅µ
c≥ s = c s −1 ⋅
ρ−ρ
si ρ = 1
K +2
1− ρ
si ρ ≠ 1
λ = λ ⋅ (1 − p K + s )
p0 =
1 s −1
∑0 c
n
+ c≥ s
n=
pn = cn ⋅ p0 pn = ρ ⋅ p n−1 qn =
L = λ ⋅W
pn
1 − pK +s
si n = 1, 2,..., s - 1 si n = s, s + 1,..., K + s si n = 0, 1,..., K + s - 1
(1 + K ⋅ ρ K +1 − ( K + 1) ⋅ ρ K ) ⋅ ρ 2 ⋅ p s −1 (1 − ρ ) 2 si ρ ≠ 1 K ⋅ ( K + 1) Lq = ⋅ p s −1 si ρ = 1 2 Lq =
5
El cálculo de W(t) para los modelos M/M/s/K, M/M/s/∞ /H y M/M/s/∞ /H con Y repuestos se explicará en el apartado 2.4.1
54
Capítulo 2:Aspectos de implementación
W = Wq +
Jorge L. Vega Valle
1
Wq =
µ
Lq λ
Cálculo de Wq(t): A = 1, S’ = 1, B = qs Desde n = s + 1 hasta K + s - 1 repetir A = A • (s • µ • t) / ( n – s ) S’ = S’ + A B = B + qn • S’ Devolver 1 – B • e –s • µ • t
2.2.2.5- Modelo M/M/1/∞ /H Intensidad = ρ = ρ ⋅ ( H − L)
ρ=
λ
c0 = 1
µ
cn = cn−1 ⋅ ρ ⋅ ( H − n + 1)
λ = λ ⋅ ( H − L)
55
si n = 1, 2,..., H
Capítulo 2:Aspectos de implementación
p0 =
Jorge L. Vega Valle
1 H
∑0 c
n
n=
L=
pn = cn ⋅ p0
si n = 1, 2,..., H
( H − n) ⋅ p n qn = ( H − L)
si n = 0, 1,..., H - 1
H
∑1 n ⋅ p
Lq = λ ⋅ W q
n
n=
W =
L
Wq = W −
λ
Cálculo de W(t):
Cálculo de Wq(t):
A = 1, S = 1, B = q0
A = 1, S = 1, B = q1
Desde n = 1 hasta (H - 1) repetir
Desde n = 2 hasta (H - 1) repetir
A = A • (µ • t) / n
A = A • (µ • t) / n
S=S+A
S=S+A
B = B + qn • S Devolver 1 – B • e
B = B + qn • S –µ • t
Devolver 1 – B • e –µ • t
56
1 µ
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.2.2.6- Modelo M/M/s/∞ /H6 Intensidad = ρ = ρ ⋅ ( H − L)
= c0 1 λ cn = cn−1 ⋅ ⋅ ( H − n + 1) si n = 1, 2,..., s n⋅µ
λ
ρ=
s⋅µ
cn = cn−1 ⋅ ρ ⋅ ( H − n + 1)
si n = s + 1,..., H
λ = λ ⋅ ( H − L)
p0 =
1 H
∑0 c
n
n=
L=
pn = cn ⋅ p0
si n = 1, 2,..., H
H n p qn = ( − ) ⋅ n ( H − L)
si n = 0, 1,..., H - 1
H
∑1 n ⋅ p
Lq = λ ⋅ W q
n
n=
W =
6
L
Wq = W −
λ
Con la restricción H ≥ s
57
1 µ
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Cálculo de Wq(t): A = 1, S’ = 1, B = qs Desde n = s + 1 hasta (H - 1) repetir A = A • (s • µ • t) / (n – s) S’ = S’ + A B = B + qn • S’ Devolver 1 – B • e –s • µ • t
2.2.2.7- Modelo M/M/s/∞ /H con Y repuestos7
Intensidad = ρ =
c0 = 1
λ
Si s ≤ Y :
s⋅µ
cn = cn−1 ⋅
cn = cn−1 ⋅ ρ ⋅ H
λ
Y+H
∑ (n − Y ) ⋅ p
n
)
n =Y
p0 =
1 Y+H
∑0 c
n
n=
p n = c n ⋅ p0
7
si n = 1, 2,..., s si n = s + 1,..., Y
cn = cn−1 ⋅ ρ ⋅ ( H + Y − n + 1) si n = Y + 1,..., Y + H Si Y ≤ s ≤ Y + H : H ⋅λ cn = cn−1 ⋅ si n = 1, 2,..., Y n⋅µ ( H + Y − n + 1) ⋅ λ cn = cn−1 ⋅ si n = Y + 1,..., s n⋅µ cn = ( H + Y − n + 1) ⋅ ρ ⋅ cn−1 si n = s + 1,..., Y + H
ρ = s⋅µ
λ = λ ⋅ (H −
H ⋅λ n⋅µ
si n = 1, 2,..., Y + H
Con la restricción H+Y ≥ s
58
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
H ⋅ pn
qn = H−
si n = 0, 1,..., Y - 1
Y +H
∑ (m − Y ) ⋅ p
m
m =Y
qn =
( H + Y − n) ⋅ pn H−
Y +H
∑ (m − Y ) ⋅ p m
si n = Y ,..., Y + H - 1
m =Y
L=
Y+H
∑1 n ⋅ p
Lq = λ ⋅ W q
n
n=
W =
L
Wq = W −
λ
1 µ
Cálculo de Wq(t): A = 1, S’ = 1, B = qs Desde n = s + 1 hasta (H - 1) repetir A = A • (s • µ • t) / (n – s) S’ = S’ + A B = B + qn • S’ Devolver 1 – B • e –s • µ • t
2.2.2.8- Modelo M/M/∞ c0 = 1
Intensidad = ρ = 0
cn = cn−1 ⋅ λ n⋅µ
59
si n = 1, 2,...
Capítulo 2:Aspectos de implementación
p0 = e
Jorge L. Vega Valle
λ − µ
p n = c n ⋅ p0
si n = 1, 2,...
L=λ µ
Lq = 0
1
Wq = 0
W =
µ
2.2.2.9- Redes de Jackson abiertas Supongamos una red de Jackson abierta con K nodos: λi = Tasa de entrada desde fuera del sistema a la cola del nodo i, i = 1, 2,..., K µi = Tasa de servicio de cada servidor del subsistema del nodo i, i = 1, 2,..., K λT
= Tasa de entrada al sistema (da igual por donde) =Σ λi
Λi
= Número efectivo de llegadas de clientes al subsistema del nodo i (da
igual por donde), i = 1, 2,..., K Matriz probabilidades: p11 p 21 M P = p i1 M p K1
p12
L
p1 j
L
p1K
p 22
L
p2 j
L
p2 K
M
O
M
O
pi 2
L
p ij
L
M
O
M
O
pK 2
L
p Kj
L
60
p iK M p KK M
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
A mayores, tenemos que:
p0 i =
λi K
∑1
i = 1, 2,..., K
pi 0 = 1 −
K
∑1 p
ij
i = 1, 2,..., K
j=
λi
i=
Denotamos con A, la siguiente matriz: 1 − p11 M A = − p1 j M −p 1K
L
− p j1
L
O
M
O
L
1 − p jj
L
O
M
O
L
− p jK
L
− p K1 M − p Kj M 1 − p KK
Si el det(A) = 0, entonces el sistema no es estacionario. Si el det(A) ≠ 0 se verifica que: Λ = inverso(A) · λ
Λ es un vector que tiene K elementos
Si Λi ≥ si · µi, para algún valor dei, entonces el sistema no es estacionario. Para que el sistema sea estacionario, debe verificarse que det(A) ≠ 0 y Λi ≤ si · µi para todos los valores de i = 1, 2,..., K. Si se verifica esta condición, podremos calcular los siguientes parámetros de salida:
61
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Obtención de parámetros relativos a un nodo:
Lq , nodo = Resolver modelo M/M/s ( Λnodo,µ nodo, s nodo ) =
c≥ s ⋅ Λnodo ⋅ p0 s nodo ⋅ µ nodo − Λnodo
(ver 2.2.2.2)
Lnodo = Lq , nodo +
Λnodo
Wnodo =
µ nodo
Lq, nodo Λnodo
+
1
Wq , nodo =
µ nodo
Lq, nodo Λnodo
Obtención de parámetros relativos a toda la red:
LT =
K
∑L
nodo
Lq , T =
nodo =1
K
∑L
WT =
q, nodo
nodo=1
LT
Wq , T =
λT
Lq , T λT
También podemos calcular fácilmente la probabilidad de que haya ni clientes en cada nodo de la red, coni = 1, 2,..., K p(n1, n2,..., nk) = p1 (n1) • p2 (n2) • … pK (nK)
siendo ni el número de clientes del nodo i.
Para obtener pi (ni), hay que operar del mismo modo que si tratase de un modelo M/M/s, con tasa de llegada Λi, tasa de servicio µi
información, ver 2.2.2.2).
62
y si servidores (para más
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.2.2.10- Redes de Jackson cerradas Supongamos una red de Jackson cerrada conK nodos y N clientes: µi = Tasa de servicio de cada servidor del subsistema del nodo i, i = 1, 2,..., K Λi = Número efectivo de llegadas de clientes al subsistema del nodo i (da igual por donde), i = 1, 2,..., K
Matriz probabilidades: (K x K elementos) p11 p 21 M P= p i1 M p K1
p12
L
p1 j
L
p1K
p 22
L
p2 j
L
p2 K
M
O
M
O
pi 2
L
p ij
L
M
O
M
O
pK 2
L
p Kj
L
p iK M p KK M
Denotamos con A, la siguiente matriz: ( (K-1) x (K-1) elementos) − p31 − p 21 1 p p 32 − − 22 A= M M − p p − 3, K −1 2, K −1
L L
− p K −1,1 − p K −1, 2
O
M
L
1 − p K −1, K −1
− p K1 − pK 2 M − p K , K −1
Denotamos con B, el siguiente vector: ( (K-1) x 1 elementos) − 1 + p11 p12 B= M p 1, K −1
63
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Λ = inverso(A) · B
Según esta fórmula Λ es un vector que tiene K-1 elementos, que se corresponden con Λ2, Λ3,..., ΛK. A mayores asignamos a Λ1 el valor de 1, con lo Λ se convierte en un vector de K elementos.
ρi =
Λi
i = 1, 2,..., K
µi
Denotamos con G, la siguiente matriz: (K x (N +1) elementos) g1 (0) g1 (1) g1 (2) g 2 (0) g 2 (1) g 2 ( 2) G = M M M g (0) g (1) g ( 2) K K K
Se verifica que: gi(0) = 1, g 1 ( n) = f 1 (n) ,
L
g1 ( N )
L
g2 (N )
O L
g K ( N ) M
para i = 1, 2,..., K para n = 0, 1,..., N
ρ in n! ρ in f i (n ) = = ai ( n) ρ in n −s si ! ⋅ si
i
si n ≤ si si n > si
Con lo que, 1 f 1 (1) f 1 ( 2) 1 g 2 (1) g 2 (2) G = M M M 1 g (1) g ( 2) K K
64
L
f1 ( N )
L
g 2 (N )
O L
g K ( N ) M
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Para calcular los valores de la matriz G que desconocemos hay que utilizar la siguiente fórmula recursiva:
g m (n) =
n
∑0 f
m
(i) ⋅ g m−1 (n − i)
m = 1, 2,..., K , n = 0, 1,..., N
i=
Ya podemos calcular fácilmente la probabilidad de que haya ni clientes en cada nodo de la red, con i = 1, 2,..., K
p n1 , n2 ,...nK =
n K ρ 1 ⋅∏ i g K ( N ) i =1 ai ( ni ) i
Cálculos relativos al último nodo:
Podemos calcular la probabilidad de que haya un determinado número de clientes (por término medio) en el último nodo:
p... i =
f K (i ) ⋅ g N ( N − i ) gK (N)
i = 0, 1,..., N
Ahora ya podemos calcular el valor deL en el último nodo:
L=
N
∑1 i ⋅ p...
i
i=
Para calcular el resto de cantidades relativas al último nodo (Lq, W y Wq), necesitamos conocer el valor de Λ K .
65
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Cálculo de Λ K :
ΛK
=0
Para i = 1, 2,..., N ≤
Si i sK ΛK
= Λ K + µK • i • p... i
ΛK
= Λ K + µK • sK • p... i
Si no Fin Si
Devolver Λ K Cálculo de W, Wq y Lq en el último nodo:
W =
L
ΛK
Wq = W −
1 µ
Lq =
L ⋅ Wq W
Para calcular estas cantidades para el resto de nodos, hay que renombrar todos los nodos de forma que el nodo de interés sea el último: Intercambiamos las posiciones en el vector µ. Intercambiamos las posiciones en el vector de servidores. Intercambiamos las posiciones en la matriz P.
66
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.3- Resolución por simulación 2.3.1- Características generales La resolución de los modelos de colas por simulación se divide en dos partes: La primera de ellas, llamada bucle de estabilización, se encarga de calcular el número de clientes que habrá en el sistema cuando entren8 en el mismo tantos como el valor del parámetro de estabilización que introduzcamos. La segunda, llamada bucle de simulación, se encarga de simular el comportamiento del sistema (teniendo en cuenta la entrada de tantos clientes como el valor del parámetro de simulación) y calcular todos los parámetros de salida del modelo. El número de clientes en el sistema al comenzar el bucle de simulación es el valor obtenido al finalizar el bucle de estabilización. A lo largo de la segunda fase de estabilización se pueden observar otras dos fases (comentadas en el apartado 1.2.5) en la obtención de un determinado parámetro de salida. Un ejemplo, tal y como se comenta en [Bos-02], lo podemos advertir en la siguiente gráfica, que corresponde a la simulación del modelo G/G/1, con distribución de llegada exponencial de parámetro λ = 2 y distribución de servicio exponencial de parámetro µ = 3.
8
En el caso de las redes cerradas, donde no hay entradas de clientes desde el exterior, tanto el parámetro de simulación como el de estabilización contabilizan las transiciones que un cliente tiene entre los nodos de la red.
67
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Figura 9- Comportamiento gráfico de la simulación
2.3.2- Implementación de los distintos modelos A continuación, presentamos en forma de pseudo-código y de modo incremental los programas de simulación (bucle de simulación) de los distintos modelos de colas. Antes de comenzar, es conveniente aclarar el significado de alguna de las variables que se usarán en los programas de simulación: c: Acumula tiempo • número de clientes que hay en el sistema. Inicialmente
vale 0. En el caso de las redes de colas será un vector con tantos elementos como el número de nodos que tengamos.
68
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
d: Acumula tiempo • número de clientes que hay en la cola. Inicialmente vale
0. En el caso de las redes de colas será un vector con tantos elementos como el número de nodos que tengamos.
2.3.2.1- Bucle de simulación del modelo G/G/1 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada, tiempo en haber una nueva salida) Si no
mínimo = tiempo que tarda en haber una nueva llegada Fin Si Si (mínimo = tiempo que tarda en haber una nueva llegada)
Aumentar el “nº de clientes simulados” en una unidad Incrementar el cronómetro en el valor mínimo Hallar el tiempo que tarda en haber una nueva llegada Incrementar el tiempo en el que hay n clientes en mínimo Si (número de clientes en el sistema = 0)
Incrementar el número de clientes en el sistema en una unidad Hallar el tiempo que tarda en haber una nueva salida Si no
Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Incrementar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo Fin Si Si no
Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema Si (número de clientes en el sistema = 0)
No puede haber una nueva salida del sistema Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema Fin Si
Disminuir el tiempo que tarda en haber una nueva llegada en mínimo Fin Si Fin Mientras
69
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.3.2.2- Bucle de simulación del modelo G/G/s9 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada, tiempo en haber una nueva salida de alguno de los servidores en los que hay un cliente) Si no
mínimo = tiempo que tarda en haber una nueva llegada Fin Si Si (mínimo = tiempo que tarda en haber una nueva llegada)
Aumentar el número de clientes simulados en una unidad Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Hallar el tiempo que tarda en haber una nueva llegada Incrementar c en mínimo • (número de clientes en el sistema) Si (número de clientes en el sistema < número de servidores)
Incrementar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Hallar el tiempo que tarda en haber una nueva salida (en uno de los servidores donde no había ningún cliente) Si no
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Incrementar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Fin Si Si no
Encontrar cual es el servidor en el cual se produce la salida de un cliente Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema)
Si número de clientes > número de servidores Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Fin Si
9
Aparecen resaltadas en color verde las diferencias respecto al modelo G/G/1
70
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Si (número de clientes en el sistema < número de servidores)
No puede haber una nueva salida del sistema del servidor en el que se ha producido la salida Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema en el servidor en el que se ha producido la salida Fin Si
Disminuir el tiempo que tarda en haber una nueva llegada en mínimo Fin Si Fin Mientras
Nota:
Como se puede observar, las diferencias entre los algoritmos de simulación para el caso de un valor de s genérico (G/G/s, G/G/s/K, G/G/s/∞ /H y G/G/s/∞ /H con Y repuestos) y el caso de un valor de s igual a 1 (G/G/1, G/G/1/K, G/G/1/∞ /H ) no son muy importantes a efectos de complejidad algorítmica. Podría haberse obviado la codificación de los casos en los que s es igual a la unidad, más se optó por su implementación para así tener un claro paralelismo entre la parte analítica y la parte de simulación. Por lo que se refiere a la parte analítica, las diferencias entre los casos en que s toma un valor genérico y en los ques es igual a 1 son muy claras, ya que es más eficiente el programar al margen este caso que el realizar los cálculos del caso genérico con el valor del número de servidores igual a la unidad. Además para algunos casos, como el M/M/1/K o M/M/1/∞ /H , su implementación proporciona una fórmula explícita para el valor de W(t), sin necesidad de emplear técnicas de integración numérica.
71
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.3.2.3- Bucle de simulación del modelo G/G/1/K10 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada, tiempo en haber una nueva salida) Si no
mínimo = tiempo que tarda en haber una nueva llegada Fin Si Si (mínimo = tiempo que tarda en haber una nueva llegada)
Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Hallar el tiempo que tarda en haber una nueva llegada Si (número de clientes en el sistema = tamaño de la c ola + 1)
Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Si no
Aumentar el “nº de clientes simulados” en una unidad Si (número de clientes en el sistema = 0)
Aumentar el número de clientes en el sistema en una unidad Hallar el tiempo que tarda en haber una nueva salida Si no
Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Aumentar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo Fin Si Fin Si Si no
Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema Si (número de clientes en el sistema = 0)
No puede haber una nueva salida del sistema
10
Aparecen resaltadas en color verde las diferencias respecto al modelo G/G/1
72
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema Fin Si
Disminuir el tiempo que tarda en haber una nueva llegada en mínimo Fin Si Fin Mientras
2.3.2.4- Bucle de simulación del modelo G/G/s/K11 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada, tiempo en haber una nueva salida de alguno de los servidores en los que hay algún cliente) Si no
mínimo = tiempo que tarda en haber una nueva llegada Fin Si Si (mínimo = tiempo que tarda en haber una nueva llegada)
Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Hallar el tiempo que tarda en haber una nueva llegada Incrementar c en mínimo • (número de clientes en el sistema) Si (número de clientes en el sistema = tamaño de cola + número de servidores)
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Si no
Aumentar el número de clientes simulados en una unidad Si (número de clientes en el sistema < número de servidores)
Incrementar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Hallar el tiempo que tarda en haber una nueva salida (en uno de los servidores donde no había ningún cliente) Si no
Incrementar d en mínimo • (nº de clientes en el sistema - nº de servidores) Incrementar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Fin Si
11
Aparecen resaltadas en color verde las diferencias respecto al modelo G/G/s
73
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Fin Si Si no
Encontrar cual es el servidor en el cual se produce la salida de un cliente Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Si número de clientes > número de servidores
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Fin Si Incrementar el tiempo en el que hay n clientes en mínimo
Disminuir en uno el número de clientes del sistema Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Si (número de clientes en el sistema < número de servidores)
No puede haber una nueva salida del servidor en el que se ha producido la salida Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema en el servidor en el que se ha producido la salida Fin Si
Disminuir el tiempo que tarda en haber una nueva llegada en mínimo Fin Si Fin Mientras
2.3.2.5- Bucle de simulación del modelo G/G/1/∞ /H Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada de algún potencial cliente, tiempo en haber una nueva salida) Si no
mínimo = tiempo que tarda en haber una nueva llegada de algún potencial cliente Fin Si Si (mínimo = tiempo que tarda en haber una nueva salida)
Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema Si (número de clientes en el sistema = 0)
No puede haber una nueva salida del sistema Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema
74
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Fin Si
Disminuir en mínimo el tiempo que tardan en llegar los potenciales clientes Hallar el tiempo que tarda en llegar un nuevo potencial cliente Si no
Aumentar el “nº de clientes simulados” en una unidad Hallar cual es el potencial cliente que tarda menos en llegar Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en mínimo el tiempo que tardan en llegar los potenciales clientes Anotar la llegada del cliente que ha tardado menos en llegar Si (número de clientes en el sistema = 0)
Incrementar el número de clientes en el sistema en una unidad Hallar el tiempo que tarda en haber una nueva salida Si no
Incrementar c en mínimo • (número de clientes en el sistema) Incrementar d en mínimo • (número de clientes en el sistema-1) Incrementar el número de clientes en el sistema en una unidad Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo Fin Si Fin Si Fin Mientras
2.3.2.6- Bucle de simulación del modelo G/G/s/∞ /H12 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada de algún potencial cliente, tiempo en haber una nueva salida de alguno de los servidores en los que hay un cliente) Si no
mínimo = tiempo que tarda en haber una nueva llegada de algún potencial cliente Fin Si Si (mínimo = tiempo que tarda en haber una nueva salida)
Encontrar cual es el servidor en el cual se produce la salida de un cliente Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema)
Si (número de clientes en el sistema > número de servidores) 12
Aparecen resaltadas en color verde las diferencias respecto al modelo G/G/1/∞ /H
75
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Fin Si Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema Disminuir en mínimo el tiempo que tarda en haber nuevas salidas de los servidores donde hay clientes
Si (número de clientes en el sistema < número de servidores) No puede haber una nueva salida del sistema de ese servidor Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema de ese servidor Fin Si
Disminuir en mínimo el tiempo que tardan en llegar los potenciales clientes Hallar el tiempo que tarda en llegar un nuevo potencial cliente Si no
Aumentar el “nº de clientes simulados” en una unidad Hallar cual es el potencial cliente que tarda menos en llegar Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Incrementar c en mínimo • (número de clientes en el sistema) Disminuir en mínimo el tiempo que tardan en llegar los potenciales clientes Anotar la llegada del cliente que ha tardado menos en llegar Disminuir el tiempo que tarda en haber una nueva salida de alguno de los servidores en los que hay un cliente en el valor mínimo Si (número de clientes en el sistema < número de servidores)
Incrementar el número de clientes en el sistema en una unidad Hallar el tiempo que tarda en haber una nueva salida de un servidor vacío Si no
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Incrementar el número de clientes en el sistema en una unidad Fin Si Fin Si Fin Mientras
76
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.3.2.7- Bucle de simulación del modelo G/G/s/∞ /H con Y repuestos13 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada de algún potencial cliente, tiempo en haber una nueva salida de alguno de los servidores en los que hay un cliente) Si no
mínimo = tiempo que tarda en haber una nueva llegada de algún potencial cliente Fin Si Si (mínimo = tiempo que tarda en haber una nueva salida)
Encontrar cual es el servidor en el cual se produce la salida de un cliente Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Si (número de clientes en el sistema > número de servidores)
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Fin Si
Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema Disminuir en mínimo el tiempo que tarda en haber nuevas salidas de los servidores donde hay clientes Si (número de clientes en el sistema < número de servidores)
No puede haber una nueva salida del sistema de ese servidor Si no
Hallar el tiempo que tarda en haber una nueva salida del sistema de ese servidor Fin Si
Disminuir en mínimo el tiempo que tardan en llegar los potenciales clientes Si (número de clientes en el sistema >= número de repuestos)
Hallar el tiempo que tarda en llegar un nuevo potencial cliente Fin Si Si no
Aumentar el “nº de clientes simulados” en una unidad Hallar cual es el potencial cliente que tarda menos en llegar Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Incrementar c en mínimo • (número de clientes en el sistema) Disminuir en mínimo el tiempo que tardan en llegar los potenciales clientes Anotar la llegada del cliente que ha tardado menos en llegar
13
Aparecen resaltadas en color verde las diferencias respecto al modelo G/G/s/∞ /H
77
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Si (número de clientes en el sistema <= número de repuestos)
Hallar el tiempo que tarda en llegar el repuesto del cliente que acaba de llegar Fin Si
Disminuir el tiempo que tarda en haber una nueva salida de alguno de los servidores en los que hay un cliente el valor mínimo Si (número de clientes en el sistema < número de servidores)
Incrementar el número de clientes en el sistema en una unidad Hallar el tiempo que tarda en haber una nueva salida de un servidor vacío Si no
Incrementar d en mínimo • (número de clientes en el sistema-número de servidores) Incrementar el número de clientes en el sistema en una unidad Fin Si Fin Si Fin Mientras
2.3.2.8- Bucle de simulación del modelo G/G/∞ Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación” Si el número de clientes en el sistema es mayor que 0
mínimo = min(tiempo en haber una nueva llegada, tiempo en haber una nueva salida de alguno de los servidores en los que hay un cliente) Si no
mínimo = tiempo que tarda en haber una nueva llegada Fin Si Si (mínimo = tiempo que tarda en haber una nueva llegada)
Aumentar el “nº de clientes simulados” en una unidad Incrementar el cronómetro en el valor mínimo Incrementar el tiempo en el que hay n clientes en mínimo Incrementar c en mínimo • (número de clientes en el sistema) Incrementar el número de clientes en el sistema en una unidad Hallar el tiempo que tarda en haber una nueva llegada Disminuir el tiempo que tarda en haber una nueva salida en el valor mínimo (en los servidores donde ya había algún cliente) Hallar el tiempo que tarda en haber una nueva salida (en un nuevo servidor donde no hay ningún cliente) Si no
Encontrar cual es el servidor en el cual se produce la salida de un cliente Incrementar el cronómetro en el valor mínimo Incrementar c en mínimo • (número de clientes en el sistema) Incrementar el tiempo en el que hay n clientes en mínimo Disminuir en uno el número de clientes del sistema
78
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Disminuir el tiempo que tarda en haber una nueva salida el valor mínimo (en los servidores donde ya había algún cliente) No puede haber una nueva salida del sistema del servidor en el que se ha producido la salida Disminuir el tiempo que tarda en haber una nueva llegada en mínimo Fin Si Fin Mientras
Para obtener el valor de los parámetros de salida14 hay que proceder como sigue:
L=
Lq =
c cronometro
d cronometro
p (n ) =
14
W =
c parametro de simulacion
Wq =
d parametro de simulacion
tiempo en el que ha habido n clientes cronometro
Referido a todos los modelos vistos hasta ahora: desde el G/G/1 hasta el G/G/∞
79
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.3.2.9 Bucle de simulación de las redes cerradas15 Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación”
Hallar cuál es el servidor (y en que nodo se encuentra) que tarda menos tiempo en atender a un determinado cliente (de aquí en adelante, ese servidor será el “servidor x”, y ese nodo será el “nodo x”) mínimo = tiempo que se tarda en atender a un cliente en el “servidor x” del “nodo x” Aumentar el “nº de clientes simulados” en una unidad Disminuir el número de clientes del “nodo x” en una unidad Disminuir en mínimo el tiempo que tardan en atender el resto de los servidores que no sean el “servidor x” a un potencial cliente. Si (nº clientes (“nodo x”) < nº servidores (“nodo x”))
No puede haber una nueva salida del “nodo x” Si no
Hallar el tiempo que tarda en producirse una nueva salida del “servidor x” del “nodo x” Fin Si
Hallar el nodo al cual ha migrado el cliente en función de la matriz de probabilidades Aumentar en uno el nº clientes (nodo al cual ha migrado) Aumentar en uno el nº entradas (nodo al cual ha migrado) Si (nº clientes (nodo al cual ha migrado) <= nº servidores (nodo al cual ha migrado))
Hallar el tiempo que tarda en producirse una nueva salida de un determinado servidor desocupado del nodo al que ha migrado Fin Si
Aumentar el cronómetro en el valor mínimo Para nodo = 1 hasta todos los nodos de la red Si (nodo = “nodo x” y nodo ~ = nodo de migración)
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo) + 1) Si (nº clientes (nodo) + 1 > nº servidores (nodo))
Aumentar el valor de d (nodo) en mínimo • (nº clientes (nodo) + 1 - nº servidores(nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo) + 1” clientes en el nodo Fin Si
15
En las redes de colas cerradas el parámetro de simulación equivale al número de transiciones de clientes entre los distintos nodos.
80
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Si (nodo ~ = “nodo x” y nodo = nodo de migración)
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo) – 1) Si (nº clientes (nodo) – 1 > nº servidores (nodo))
Aumentar el valor de d (nodo) en mínimo • (nº clientes (nodo) – 1 - nº servidores(nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo) – 1” clientes en el nodo
Si no
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo)) Si (nº clientes (nodo) > nº servidores (nodo))
Aumentar el valor de d (nodo) en mínimo • (nº clientes (nodo) – nº servidores(nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo)” clientes en el nodo Fin Si Fin Para Fin Mientras
Inicialmente, los clientes se distribuyen de la siguiente forma:
n º clientes red ⋅ (nº servidores (nodo))-1 Nodos de la red cerrada -1 (nº servidore s (i)) ∑ i =1
n º clientes ( nodo) =
Según esta fórmula es posible que quede algún cliente sin asignar, es decir, que tras esta distribución inicial, la suma del número de clientes que hay en todos los nodos sea inferior al número de clientes que ha introducido el usuario. Para que esto no ocurra, se inicia un reparto secuencial de los “clientes residuales 16 ”, empezando por el primer nodo.
Nótese que las actualizaciones de los valores de c, d y el tiempo que un cliente está en un determinado nodo son “a posteriori”. Así, si tenemos la siguiente situación:
16
Se corresponden con el número de clientes que ha introducido el usuario menos la suma del número de clientes en los nodos de la red tras aplicar la fórmula inicial
81
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Nodo de migración
Nodo x
3 clientes antes de migrar
5 clientes antes de migrar
2 clientes despues de migrar
6 clientes despues de migrar
Otros Nodos
Figura 10- Actualización de los valores de " c" y "d" en las redes cerradas
Al entrar en el bucle “para” que recorre todos los nodos haremos las siguientes operaciones (en los nodos etiquetados como “nodo x” y “nodo de migración” de la anterior figura): Nodo x c (nodo x) = c (nodo x) + mínimo • 3
Si ( 3 > nº de servidores (nodo x) ) d (nodo x) = d (nodo x) + mínimo • ( 3 – nº servidores (nodo x) )
Tiempo (3 clientes, nodo x) = Tiempo (3 clientes, nodo x) + mínimo Nodo de migración c (nodo migración) = c (nodo migración) + mínimo • 5
Si ( 5 > nº de servidores (nodo migración) ) d
(nodo
migración)
=
d
(nodo
migración)
+
mínimo
•
( 5 – nº servidores (nodo migración) ) Tiempo (5 clientes, nodo migración) = Tiempo (5 clientes, nodo migración) + mínimo
82
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Para obtener el valor de los parámetros de salida hay que proceder como sigue:
c ( nodo) cronometro
L (nodo) =
Lq (nodo) =
d (nodo) cronometro
Lq , T =
Nodos de la red cerrada
∑ i 1 =
W (nodo) =
c (nodo) nº de entradas (nodo)
Wq ( nodo) =
d (nodo) n º de entradas (nodo)
d (i) cronometro
p ( n, nodo) = tiempo en el que hay n clientes en un determinado nodo cronometro
2.3.2.10- Bucle de simulación de las redes abiertas Mientras “nº de clientes simulados” < “nº de clientes totales en la simulación”
mínimo = min(tiempo que tarda en salir un determinado cliente de un “servidor x” de un “nodo x”, tiempo que tarda en llegar un potencial cliente a un “nodo y”) Si (mínimo = tiempo de llegada de un potencial cliente)
Aumentar en una unidad el “nº de clientes simulados” Aumentar en una unidad el nº de entradas al “nodo y” Aumentar en una unidad el nº de clientes del “nodo y” Disminuir en mínimo el tiempo de llegada de potenciales clientes Hallar el tiempo de llegada de un nuevo potencial cliente al “nodo y” Disminuir en mínimo el tiempo de salida de clientes en los servidores de aquellos nodos donde exista un cliente
83
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Si (nº clientes (“nodo y”) < nº servidores (“nodo y”))
Generar el tiempo de salida de un servidor del “nodo y” que estaba desocupado Fin Si
Aumentar la variable cronómetro en el valor mínimo Para nodo = 1 hasta t odos los nodos de la red
Si (nodo = “nodo y”) Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo) – 1) Si (nº clientes (nodo) – 1 > nº servidores (nodo))
Aumentar d (nodo) en mínimo • (nº clientes (nodo) – 1 – nº servidores (nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo) – 1” clientes en el nodo Si no
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo)) Si (nº clientes (nodo) > nº servidores (nodo))
Aumentar d (nodo) en mínimo • (nº clientes (nodo) – nº servidores (nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo)” clientes en el nodo Fin Si Fin Para Si no
Disminuir en uno el número de clientes del “nodo x” Disminuir en mínimo el tiempo de salida de clientes en los servidores de aquellos nodos donde exista un cliente Si (nº clientes (“nodo x”) < nº servidores (“nodo x”))
No puede haber una nueva salida del “servidor x” del “nodo x” Si no
Calcular el tiempo que tarda en producirse una nueva salida en el “servidor x” del “nodo x” Fin Si
Disminuir en mínimo el tiempo de llegada de potenciales clientes Hallar el nodo (o exterior) al cual migra el cliente en función de la matriz de probabilidades Si (nodo de migración ~= exterior del sistema)
Aumentar en uno el nº clientes (nodo al cual ha migrado) Aumentar en uno el nº entradas (nodo al cual ha migrado) Si (nº clientes (nodo de migración) <= nº servidores (nodo de migración)) Hallar el tiempo de salida de este nuevo cliente que ocupa un servidor que estaba desocupado en el nodo de migración Fin Si Fin Si
84
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
Aumentar el cronómetro en el valor mínimo Para nodo = 1 hasta t odos los nodos de la red Si (nodo = “nodo x” y nodo ~ = nodo de migración)
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo) + 1) Si (nº clientes (nodo) + 1 > nº servidores (nodo))
Aumentar el valor de d (nodo) en mínimo • (nº clientes (nodo) + 1 - nº servidores(nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo) + 1” clientes en nodo Fin Si Si (nodo ~ = “nodo x” y nodo = nodo de migración)
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo) – 1) Si (nº clientes (nodo) – 1 > nº servidores (nodo))
Aumentar el valor de d (nodo) en mínimo • (nº clientes (nodo) – 1 - nº servidores(nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo) – 1” clientes en nodo Si no
Aumentar el valor de c (nodo) en mínimo • (nº clientes (nodo)) Si (nº clientes (nodo) > nº servidores (nodo))
Aumentar el valor de d (nodo) en mínimo • (nº clientes (nodo) – nº servidores(nodo)) Fin Si
Aumentar en mínimo el tiempo en el que hay “nº clientes (nodo)” clientes en nodo Fin Si Fin Para Fin Mientras
Para obtener el valor de los parámetros de salida hay que proceder como sigue:
L (nodo) =
c ( nodo) cronometro
Lq (nodo) =
d (nodo) cronometro
85
W (nodo) =
c (nodo) nº de entradas (nodo)
Wq (nodo) =
d (nodo) nº de entradas (nodo)
Capítulo 2:Aspectos de implementación
Nodos de la red cerrada
LT =
∑ i =1
L q, T =
Jorge L. Vega Valle
c (i) cronometro
Nodos de la red cerrada
∑ i =1
p ( n, nodo) =
WT =
d (i) cronometro
W q, T =
LT
nº de clientes simulacion cronometro
L q, T
nº de clientes simulacion cronometro
tiempo en el que hay n clientes en un determinado nodo cronometro
86
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
2.4- Detalles de la implementación En este apartado se explica con mayor detalle alguna de las características de la aplicación, así como algunos de los problemas surgidos a lo largo de la implementación y la forma de solucionarlos.
2.4.1- Cálculo de W(t) en los modelos en los que no hay una solución analítica Tanto para el modelo M/M/s/K (ver apartado 2.2.2.4), como para los modelos M/M/s/∞ /H (apartado 2.2.2.6) y M/M/s/∞ /H con Y repuestos (apartado 2.2.2.7) no
se ha proporcionado una fórmula explícita para el cálculo de W(t). Esto es debido a que dicha fórmula es excesivamente compleja, y la mayoría de los libros ni siquiera la consideran. Por todo ello, hemos tenido que optar por otra vía para calcular el valor de W(t). En todos los modelos se verifica que:
W (t ) =
t
∫0 W
q
(t − s ) ⋅ µ ⋅ e − µ ⋅ s ds
La resolución numérica exacta de dicha integral es bastante complicada, por lo que, para hallar el valor de W(t) en los modelos en los que no tenemos una fórmula exacta, hemos utilizado la implementación que MATLAB realiza del método de Simpson17 (función quad) de integración numérica. Pese a lo que se pueda pensar en un primer momento, el hecho de utilizar un método de integración numérica para hallar el valor de W(t) apenas tiene influencia en el tiempo de ejecución del algoritmo y el sistema nos proporciona una respuesta 17
Puede observarse con detalle en [Gar-98]
87
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
de forma muy rápida. Incluso para el caso de tener que mostrar la evolución de la función W(t) en un rango de valores (esto se hace resolviendo la integral para cien valores equiespaciados en dicho rango, es decir realizando cien aproximaciones mediante integración numérica) se resuelve muy rápidamente.
2.4.2- Problemas de desbordamiento en el cálculo Debido a como se resuelven analíticamente muchos de los modelos, se pueden plantear problemas de desbordamiento en el cálculo de algunos parámetros. Un ejemplo puede ser el siguiente: Modelo M/M/1/K con K = 300, λ = 100 y µ = 1 En este caso, ρ = 100. Recordemos la fórmula de L, para el modelo M/M/1/K:
L=
L=
ρ
1− ρ
−
( K + 2) ⋅ ρ K + 2 1 − ρ K +2
, sustituyendo tenemos que
100 302 ⋅ 100 302 − , lo que provoca un desbordamiento en el cálculo, ya que − 99 1 − 100 302
100302 es un número superior al mayor número real que puede manejar MATLAB (1.7976 • 10308) Este problema puede darse en otros modelos en función de los valores de los parámetros de entrada que introduzca el usuario. Para solucionar el problema, le indicamos al usuario por medio de un mensaje de error que se ha producido un
88
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
desbordamiento en la realización de los cálculos y debe cambiar los valores de los parámetros de entrada.
2.4.3- Problemas de precisión en el cálculo A lo largo del desarrollo de la aplicación han aparecido varios problemas relativos a la precisión en el cálculo de las distintas variables. A continuación se citan los más importantes y la manera de resolverlos: En las redes cerradas se debe verificar la condición de que la suma de las probabilidades de transición de un nodo concreto a todos los nodos de la red debe ser igual a 1. Pues bien, si por ejemplo en MATLAB hacemos 0.4 + 0.3 + 0.2 + 0.1, esto da como resultado un valor distinto de 1!!! Para resolver este problema hemos sustituido la condición de que la suma sea igual a 1, por la condición de que la suma se encuentre entre 0.999999 y 1.000001. En muchos modelos, el cálculo de Wq y Lq se realiza según las siguientes fórmulas: Wq = W −
1
L q = λ ⋅ Wq
µ
Esto puede provocar que debido a las múltiples operaciones que hay que realizar para el cálculo de L (normalmente son varias sumas y productos), se puede obtener en algunos casos un valor para W ligeramente inferior a 1/µ, lo que provocaría que tanto el valor de Wq, como de Lq fuesen negativos, en vez de ser igual a cero. Para solucionar este problema, asignamos directamente el valor cero a dichas variables en caso de que se de esta situación. Debido a que, para el cálculo analítico de W(t) en los modelos M/M/s/K, M/M/s/∞ /H y M/M/s/∞ /H con Y repuestos se utilizan técnicas de integración
89
Capítulo 2:Aspectos de implementación
Jorge L. Vega Valle
numérica, puede producirse que el valor de W(t) sea, en algunos casos, ligeramente superior a 1. Para solucionar este problema se asigna directamente un valor igual a 1 en caso de que ocurra. Algunas veces al usar las técnicas de integración numérica pueden aparecer los siguientes mensajes de aviso en la línea de comandos de MATLAB :
‘Warning:
Minimum step size reached; singularity possible ’
Este mensaje se produce cuando al intentar evaluar el valor de la integral y debido a las características de la misma, MATLAB no puede garantizar que el error cometido en dicha evaluación sea inferior al requerido, que en nuestro caso es de una milésima. ‘Warning:
Infinite or Not-a-Number function value
encountered’
Este mensaje indica que se ha producido un error de overflow (desbordamiento de cálculo) o una división por cero durante la evaluación de la integral. Para evitar estos mensajes de aviso basta con escribir lo siguiente (esto sólo es válido con la versión 6.5 de MATLAB®) en la línea de comandos: ‘warning
off MATLAB:quad:MinStepSize ’
‘warning
off MATLAB:quad:ImproperFcnValue ’
90
ó
Capítulo 3: Evaluación de la aplicación Tolle numerorum omnibus rebus et omnia pereunt “Quita el número de las cosas y todas se destruirán” (San Isidoro de Sevilla, obispo y teólogo de la España visigoda, siglo VIII)
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
3.1- Validación de los resultados de la simulación Para comprobar la corrección de los resultados obtenidos, compararemos los resultados obtenidos mediante simulación con los de la solución analítica correspondiente al mismo modelo. Nótese que para poder validar los resultados de la simulación con los de la solución analítica deberemos de considerar tanto la distribución de llegada como la de servicio de tipo exponencial. La validación de los resultados de la simulación se enfocará estudiando la influencia de los siguientes parámetros: Parámetro de estabilización y número de clientes de la simulación Intensidad de tráfico Modelo de colas
3.1.1- Variación del parámetro de estabilización y del número de clientes de la simulación Para observar la influencia del parámetro de estabilización y del número de clientes de la simulación en el resultado del algoritmo de simulación se han analizado los dos siguientes casos: Primer caso
Modelo de colas: Simulación del modelo G/G/s/K Distribución de llegada: Exponencial, con λ = 10 Distribución de servicio: Exponencial, con µ = 5 Número de servidores: s = 3
92
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Tamaño de la cola: K = 3
Nº Clientes Simulación
. m á r a P
. b a ts E
100
1000
5000
10000
50000
100000
8.582 %
7.017 %
1.415 %
1.643 %
0.481 %
0.463 %
1000 19.235 % 3.998 %
3.520 %
0.463 %
0.676 %
0.106 %
10000 11.714 % 5.213 %
2.281 %
0.885 %
0.336 %
0.339 %
10000
50000
100000
3.935 %
0.680 %
0.338 %
100
Segundo caso
Modelo de colas: Simulación del modelo G/G/s/K Distribución de llegada: Exponencial, con λ = 200 Distribución de servicio: Exponencial, con µ = 5 Número de servidores: s = 30 Tamaño de la cola: K = 300
Nº Clientes Simulación 100 . m á r a P
. b a ts E
100
1000
5000
83.290 % 52.723 % 6.986 %
1000 17.209 % 3.342 %
0.257 %
0.237 %
0.040 %
0.040 %
10000 0.545 %
0.030 %
0.043 %
0.021 %
0.025 %
0.130 %
En las tablas anteriores, se muestra el error medio (en tanto por ciento) de cinco muestras en el cálculo de L.
93
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Como podemos observar, los resultados del primer caso son muy diferentes a los del segundo. Así, mientras en el primer caso apenas tiene influencia el aumento del valor del parámetro de estabilización, para el segundo caso podemos ver como un aumento del parámetro de estabilización tiene una gran influencia en el error obtenido. Esto es debido a que en el primer caso no es muy importante el estado del cual se parta para el cálculo de L, mientras que en el segundo caso es crucial. En ambos casos podemos ver como cuanto mayor es el número de clientes de la simulación, menor es el error que se comete. La disminución del error es muy pronunciada al principio y más lenta a medida que ya estamos muy cerca de la solución óptima (nótese que el pasar de 50.000 a 100.000 clientes en la simulación apenas tiene influencia en el valor del error). Recomendación práctica
Debido a que el parámetro de estabilización puede tener una gran influencia en el modelo (dependiendo de las características del mismo) y a que un valor muy elevado para el número de clientes de la simulación no hace sino ralentizar el tiempo de respuesta, los siguientes valores pueden considerarse como los más adecuados para obtener lo más rápidamente posible una respuesta con un pequeño valor para el error: Parámetro de estabilización = 5.000 Número de clientes de la simulación = 10.000
94
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
3.1.2- Variación de la intensidad de tráfico Para estudiar la influencia de la intensidad de tráfico en el error cometido por la simulación vamos a analizar tres posibles situaciones: Primer caso
Modelo de colas: Simulación del modelo G/G/s Distribución de llegada: Exponencial Distribución de servicio: Exponencial, con µ = 5 Número de servidores: s = 2 Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000
λent
Intensidad
Error
1
0.1
0.983 %
5
0.5
1.868 %
7
0.7
3.771 %
9
0.9
11.326 %
9.5
0.95
25.2192 %
Segundo caso
Modelo de colas: Simulación del modelo G/G/s/∞ /H Distribución de llegada: Exponencial Distribución de servicio: Exponencial, con µ = 5 Número de servidores: s = 2 Población potencial: H = 5
95
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000
λent
Intensidad
Error
0.2
0.096137
2.001 %
1.25
0.48824
0.943 %
2.25
0.71267
1.662 %
4
0.89162
0.673 %
6
0.9588
0.750 %
Tercer caso
Modelo de colas: Simulación del modelo G/G/s/∞ /H Distribución de llegada: Exponencial Distribución de servicio: Exponencial, con µ = 5 Número de servidores: s = 2 Población potencial: H = 500 Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000
λent
Intensidad
Error
0.002
0.09996
1.653 %
0.01003
0.50017
3.357 %
0.01408 0.01835
0.70018 0.90181
4.476 % 7.201 %
0.01977
0.95719
17.8428 %
96
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Los resultados obtenidos en las tablas anteriores muestran el error medio cometido (por cinco muestras) en el cálculo de L. Como podemos ver en las tablas anteriores, dependiendo de las características del modelo que tengamos, el valor de la intensidad de tráfico puede tener una gran influencia en el error cometido. Para el primer y tercer caso, vemos como al aumentar el valor de la intensidad de tráfico mayor es el error cometido y, a medida que ésta se acerca más a uno, el error aumenta mucho más. Esto es debido a que en estos casos, cuanto mayor sea el valor de la intensidad de tráfico, mayor será también la dependencia que existe entre los datos de la muestra obtenida por simulación. Por otra parte, para el segundo caso no ocurre lo mismo, ya que, al tener un valor de H tan bajo, la dependencia entre las observaciones se atenúa (el tiempo de estancia en el sistema de un cliente y del siguiente a él no pueden tener una correlación muy alta)
3.1.3- Variación del modelo de colas En este apartado estudiaremos la influencia que tiene la elección del modelo de colas en el error cometido. Debido a que, como hemos observado en el apartado anterior, la intensidad de tráfico tiene influencia en el valor del error (al menos, para algunos casos), consideraremos modelos cuya intensidad de tráfico sea razonablemente baja. Para observar la influencia del parámetro modelo de colas en el resultado del algoritmo de simulación se han fijado los siguientes parámetros: Distribución de llegada: Exponencial Distribución de servicio: Exponencial
97
Capítulo 3: Evaluación de la aplicación
Intensidad de tráfico
18
Jorge L. Vega Valle
: Baja
Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000 Los resultados obtenidos pueden verse en la siguiente tabla: (el error que se muestra corresponde a la media de cinco muestras en el cálculo de L) λent
µsal
Modelo de colas
Parámetros del
Intensidad
Error
1
10
G/G/1
0.1
0.875 %
1
5
G/G/s
= 2s
0.1
1.099 %
2.75
25
G/G/1/K
=K5
0.11
1.209 %
2
5
G/G/s/K
s = 4, =K5
0.1
1.529 %
2.1
40
G/G/1/∞ /H
=H2
0.099515
0.624 %
2
28
G/G/s/∞ /H
s = 4,=H6
0.1
0.963 %
2
40
2
5
modelo
G/G/s/∞ /H con Y
repuestos
s = 4, H = 8, Y= 2
0.099893 0
G/G/∞
0.806 % 1.307 %
K = 2,
Redes Jackson abiertas
λ1 = 2, µ1 =
5, s1 = 5,
p11 = 0.2, p12 = 0.1, λ2 =
1.039 %
3, µ2 = 4, s2 = 10,
p21 = 0.1, p22 = 0.1 K = 2, N = 4,
Redes Jackson cerradas
µ1 = 2, s1 = 5, p11 = 0.1, p12 = 0.9, µ2 = 3, s2 = 5, p21 = 0.4, p22 = 0.6
18
No definida como tal para los modelos de redes de colas
98
1.126 %
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Como se puede observar en la tabla anterior, los resultados de la simulación son muy parecidos a los del modelo analítico correspondiente. De hecho, un estudio más completo muestra que solamente no es así cuando nos encontramos en modelos con una intensidad de tráfico elevada en donde haya una gran dependencia entre los datos de la muestra obtenida por simulación.
99
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
3.2- Validación de la simulación de las distribuciones de probabilidad19 Para verificar el correcto funcionamiento de la simulación de las distintas distribuciones (ver [Cao-98] o [Dev-86]) que se han empleado en la aplicación, se ha usado la fórmula de Pollaczek-Khintchine, explicada con detalle en [Cao-02], que nos permite obtener el valor de L en modelos de la forma M/G/1.
L=ρ+
ρ 2 + λ2 ⋅ σ S2
2 ⋅ (1 − ρ )
donde: • ρ es la intensidad de tráfico, y en este caso equivale a:
ρ=
λ µ
=
λ
= λ ⋅ Media ( Distrib de Servicio )
1 Media ( Distribución de Servicio)
Además debe cumplirse queρ < 1, para así garantizar la estacionariedad. • λ es el número medio de llegadas de clientes al sistema por unidad de tiempo
(recordemos que para poder aplicar la fórmula de Pollaczek-Khintchine la distribución de llegada de clientes al sistema debe ser siempre exponencial). • σ2S es la varianza de la distribución de servicio. 19
Ver Apéndice A para observar la definición de dichas distribuciones
100
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Los resultados de las distintas pruebas son los siguientes:20 Distribución exponencial
El valor de ρ tiene una gran influencia en los resultados obtenidos durante la simulación. Así, el error cometido es mucho mayor cuando ρ toma valores cercanos a uno que cuando se encuentra alejado de este valor. Como ejemplo, vamos a comentar los resultados obtenidos para el caso de la distribución exponencial: Caso 1: Intensidad de tráfico baja ( = 0.15) λ=3 µ=
20
ρ= 0.15 σ
2
S
= 0.0025 L(Pollaczek-Khintchine) = 0.17647
L (Simulación)
Error cometido (en %)
Prueba nº1
0.17976
1.8643
Prueba nº2
0.17678
0.1757
Prueba nº3
0.18051
2.2893
Prueba nº4
0.17397
1.4167
Prueba nº5
0.17623
0.1360
Error Medio
1.1764
Media(L(Simulación))
Error de la Media (en %)
0.17745
0.5553
20
Todas las pruebas se han realizado con un valor igual a 1.000 en parámetro de estabilización y con un valor igual a 10.000 en número de clientes de la simulación
101
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Caso 2: Intensidad de tráfico media ( = 0.60) λ=3 µ=5 ρ= σ
2
S
0.6 = 0.04 L(Pollaczek-Khintchine) = 1.5
(Simulación) L
Error cometido (en %)
Prueba nº1
1.4866
0.8933
Prueba nº2
1.4456
3.6266
Prueba nº3
1.4892
0.7200
Prueba nº4
1.6302
8.6800
Prueba nº5
1.4776
1.4933
Error Medio
3.0826
Media(L(Simulación))
Error de la Media (en %)
1.50584
0.3893
Caso 3: Intensidad de tráfico alta ( = 0.90) λ=3 µ=
10/3
ρ=
0.9
σ
2
S
= 0.09
102
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
L(Pollaczek-Khintchine) = 9
L (Simulación)
Error cometido (en %)
Prueba nº1
10.6167
17.9634
Prueba nº2
8.5961
4.4878
Prueba nº3
9.3125
3.4723
Prueba nº4
8.4652
5.9423
Prueba nº5
12.0142
33.4912
Error Medio
13.0714
Media(L(Simulación))
Error de la Media (en %)
9.8009
8.8989
En el siguiente gráfico se puede observar de forma clara lo comentado anteriormente. Cuando ρ toma valores elevados el error medio es mucho mayor (13%) que cuando toma valores medios (3%) o bajos (1%).
Influencia de la intensidad de tráfico en el error cometido Error (%)
35,0000 30,0000 25,0000 20,0000
Ro = 0.15 (bajo)
15,0000
Ro = 0.6 (medio) Ro = 0.9 (alto)
10,0000 5,0000 0,0000
Prueba Prueba Prueba Prueba Prueba 1 2 3 4 5
103
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución uniforme λ=1
a = 0.2 b = 0.8 ρ = 0.5 σ
2
S
= 0.03 L(Pollaczek-Khintchine) = 0.78
(Simulación) L
Error cometido (en %)
Prueba nº1
0.79093
1.4013
Prueba nº2
0.78149
0.1910
Prueba nº3
0.77515
0.6218
Prueba nº4
0.75774
2.8538
Prueba nº5
0.78549
0.7038
Error Medio
1.1543
Media(L(Simulación))
Error de la Media (en %)
0.77816
0.2359
104
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución determinista λ=2
d = 0.3 ρ = 0.6 σ
2
S
=0 L(Pollaczek-Khintchine) = 1.05
(Simulación) L
Error cometido (en %)
Prueba nº1
1.0164
3.2000
Prueba nº2
1.0245
2.4286
Prueba nº3
1.1154
6.2286
Prueba nº4
1.0288
2.0190
Prueba nº5
1.0466
0.3238
Error Medio
2.8400
Media(L(Simulación))
Error de la Media (en %)
1.04634
0.3486
105
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución gamma λ
= 0.5
p=2 a=3 ρ = 1/3 σ
2
S
= 2/9
≃
0.2222
L(Pollaczek-Khintchine) = 11/24
(Simulación) L
≃
0.4583333
Error cometido (en %)
Prueba nº1
0.4621
0.8218
Prueba nº2
0.4751
3.6582
Prueba nº3
0.4529
1.1855
Prueba nº4 Prueba nº5
0.4620 0.4528
0.8000 1.2073
Error Medio
1.5346
Media(L(Simulación))
Error de la Media (en %)
0.4610
0.5818
106
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución beta λ=2
p=1 q=5 k=1 ρ = 1/3 σ
2
S
= 5/252
L(Pollaczek-Khintchine)= 120/252
(Simulación) L
≃
0.47619
Error cometido (en %)
Prueba nº1
0.4795
0.6950
Prueba nº2
0.4759
0.0610
Prueba nº3 Prueba nº4
0.4836 0.4872
1.5560 2.3120
Prueba nº5
0.4824
1.3040
Error Medio
1.1856
Media(L(Simulación))
Error de la Media (en %)
0.4817
1.1570
107
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución lognormal λ
= 0.2
µ=1 σ = 0.2 ρ = 0.2 * e σ
2
S
1.02
≃
= e2.08 - e2.04
0.5546
≃
0.3139
L(Pollaczek-Khintchine)= 0.2 * e1.02 + (0.02 * e2.08) / (1 – 0.2 * e1.02)
(Simulación) L
≃
0.9141
Error cometido (en %)
Prueba nº1
0.9212
0.7769
Prueba nº2
0.9038
1.1267
Prueba nº3
0.8895
2.6910
Prueba nº4
0.8789
3.8506
Prueba nº5
0.9090
0.5578
Error Medio
1.8006
Media(L(Simulación))
Error de la Media (en %)
0.9005
1.4877
108
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución normal λ
= 0.2
µ=4 σ = 0.5 ρ = 0.8 σ
2
S
= 0.25 L(Pollaczek-Khintchine) = 2.425
L (Simulación)
Error cometido (en %)
Prueba nº1
2.3993
1.0598
Prueba nº2
2.6071
7.5093
Prueba nº3
2.3044
4.9732
Prueba nº4
2.5243
4.0948
Prueba nº5
2.3460
3.2577
Error Medio
4.1790
Media(L(Simulación))
Error de la Media (en %)
2.4362
0.4627
109
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Distribución de Weibull λllegada
=1
λWeibull = 10 α=
0.5
ρ=
0.2
σ
2
S
= 0.2 L(Pollaczek-Khintchine) = 0.35
L (Simulación)
Error cometido (En %)
Prueba nº1
0.3703
5.8000
Prueba nº2
0.3530
0.8571
Prueba nº3
0.3546
1.3143
Prueba nº4
0.3595
2.7143
Prueba nº5
0.3657
4.4857
Error Medio
3.0343
Media(L(Simulación))
Error de la Media (En %)
0.3606
3.0286
Como se puede observar en los resultados de todas las distribuciones, el error cometido es siempre inferior a un 5%, por lo que la simulación de las distintas funciones de distribución funciona correctamente. Solamente hay que tener precaución con los resultados de la simulación cuando los valores de la intensidad de tráfico son cercanos a uno.
110
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
3.3- Influencia de la variación de algunos parámetros en el tiempo de ejecución de la simulación En este apartado se analizará la influencia que tiene la variación de algunos parámetros en el tiempo de ejecución de los algoritmos de simulación. En el caso de la resolución analítica la respuesta es inmediata en casi todos los casos (salvo en las redes de colas cerradas que posean un número de nodos y clientes elevados). Todas las pruebas han sido realizadas en el siguiente equipo: Pentium III a 933 Mhz 384 Mb de memoria RAM Carga mínima del sistema (máquina aislada).
Los factores a considerar son los siguientes: Distribución de probabilidad Intensidad de tráfico Parámetro de estabilización y número de clientes de la simulación Modelo de colas
3.3.1- Variación de la distribución de probabilidad Para observar la influencia de una determinada distribución de probabilidad en el tiempo de ejecución del algoritmo de simulación se han fijado los siguientes parámetros:
111
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Modelo de colas: Simulación del modelo G/G/1 Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000 Distribución de llegada: Exponencial La distribución de servicio varía en función de que distribución de probabilidad estemos estudiando. Los resultados obtenidos pueden verse en la siguiente tabla21: λent
Distrib. Servicio
Parámetros de la
Intensidad
Tiempo
distribución
3
EXPONENCIAL
= 0.6µ
0.6
3.6032
1
UNIFORME
a = 0.2, = 0.8 b
0.5
5.8804
2 2
DETERMINISTA GAMMA
= 0.3d a = =5,1p
0.6 0.4
3.5650 7.3848
2
BETA
p = 1, q = 5,= k1
0.2
LOGNORMAL
µ = 1, = 0.2 σ
0.3333
17.7636
0.5546
5.8624
0.2
NORMAL
µ ==4,0.5 σ
0.8
6.1750
1
WEIBULL
λ = 10, = 0.5 α
0.2
6.2448
Como conclusión, se puede observar que la distribución que se elija tiene una gran influencia en el tiempo de ejecución del algoritmo; por ejemplo, el tiempo de ejecución del algoritmo con la distribución beta es casi cinco veces superior al de la distribución determinista. Esto es debido a que la generación de números de una determinada distribución (ver 2.1.4) es más compleja en unos casos que en otros.
21
Puede observarse que todos los modelos son estacionarios ( ρ < 1). Los tiempos de ejecución que se muestran en todas las tablas son en segundos.
112
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
3.3.2- Variación de la intensidad de tráfico Para observar la influencia de la intensidad de tráfico en el tiempo de ejecución del algoritmo de simulación se han fijado los siguientes parámetros: Modelo de colas: Simulación del modelo G/G/1 Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000 Distribución de llegada: Exponencial Distribución de servicio: Exponencial, con µ = 10 El valor de λ varía en función del valor de la intensidad de tráfico que estemos estudiando. Los resultados obtenidos pueden verse en la siguiente tabla: λent
Intensidad
Tiempo
0.1
0.01
3.5230
1
0.1
3.5570
3
0.3
3.6110
5
0.5
3.6736
7
0.7
3.7234
9
0.9
3.7654
9.5
0.95
3.7854
9.9
0.99
3.8056
Como conclusión, se puede observar que el valor de la intensidad de tráfico tiene una ligera influencia en el tiempo de ejecución del algoritmo. Esto es debido a que para la finalización del bucle de simulación (y para la finalización del bucle de
113
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
estabilización), se requiere que el número de llegadas, que inicialmente vale cero, sea igual al número de clientes de la simulación (en el caso del bucle de estabilización, igual al parámetro de estabilización). El número de llegadas aumenta cuando no hay ningún cliente en el sistema (debe de producirse de forma forzosa una llegada) o cuando hay algún cliente y se produce una llegada. La intensidad de tráfico tiene una gran influencia en ambos casos, pues, el primero de ellos se produce con mayor frecuencia cuanto menor es su valor, mientras que el segundo de los casos se produce con mayor frecuencia cuanto mayor es su valor.
3.3.3- Variación del parámetro de estabilización y número de clientes de la simulación Para observar la influencia del parámetro de estabilización y del número de clientes de la simulación en el tiempo de ejecución del algoritmo de simulación se han fijado los siguientes parámetros: Modelo de colas: Simulación del modelo G/G/1 Distribución de llegada: Exponencial, con λ = 3 Distribución de servicio: Exponencial, con µ = 5 Los resultados obtenidos pueden verse en la siguiente tabla:
114
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Nº Clientes Simulación
e d o r t
n ó i c a
e m á r a P
ilz b a ts E
1000
5000
10000
20000
50000
100000
10
1.913
2.614
3.555
5.417
10.966
20.159
100
1.893
2.634
3.565
5.408
10.986
20.169
500
1.953
2.704
3.626
5.518
11.086
20.229
1000
2.003
2.764
3.686
5.558
11.186
20.329
5000
2.544
3.355
4.226
6.138
11.617
20.850
10000
3.304
3.976
4.937
6.760
12.418
21.501
Como conclusiones podemos citar las siguientes: El aumento del parámetro de estabilización tiene especial influencia en el tiempo de ejecución del algoritmo de simulación cuando el número de clientes de la simulación es bajo. Podemos observar, que con 1.000 clientes de simulación, la duración del algoritmo en función del parámetro de estabilización fluctúa un 72%, con 10.000 clientes un 39%, y con 100.000 clientes sólo oscila un 7%. La complejidad del algoritmo de simulación es:
Ο
(n). La comprobación
puede verse en la siguiente tabla, realizada para el caso en el cual el parámetro de estabilización vale 10.000 (en otro caso ocurriría algo muy similar): Cociente entre el t. ejecución y las siguientes funciones
s e t n e il C º N
n ió c a l u m i S
log (n)
N
n · log(n)
n2
1000
0.826
3.304·10-3
1.101·10-3
3.304·10-6
5000
1.075
7.952·10-4
2.150·10-4
1.590·10-7
-4
-4
-8
10000 20000
1.234 1.572
4.937·10 3.380·10-4
1.23410 7.859·10-5
4.937·10 1.690·10-8
50000
2.643
2.484·10-4
5.285·10-5
4.967·10-9
100000
4.300
2.150·10-4
4.300·10-5
2.150·10-9
115
Capítulo 3: Evaluación de la aplicación
Jorge L. Vega Valle
Para el caso de la función log(n), la sucesión parece ser divergente. Se trata de una función “subestimada22 ”. Para n2 podemos ver que la sucesión aparenta converger a cero. Se trata de un función “sobrestimada”. Tanto al considerar la función n, como la función n·log(n), las sucesiones parecer converger a un valor distinto de cero, por lo que se trata de funciones “ajustadas”, aunque la primera con una mayor precisión que la segunda.
3.3.4- Variación del modelo de colas Para observar la influencia del modelo de colas en el tiempo de ejecución del algoritmo de simulación se han fijado los siguientes parámetros: Parámetro de estabilización: 1.000 Número de clientes de la simulación: 10.000 Distribución de llegada: Exponencial Distribución de servicio: Exponencial Intensidad de tráfico23: Igual ó muy cercana a 0.4 Los resultados obtenidos pueden verse en la siguiente tabla:
22 23
Estos términos se definen con detalle en [Wei-95] Salvo en el modelo G/G/∞ en donde la intensidad de tráfico es siempre igual a cero
116
Capítulo 3: Evaluación de la aplicación
λent
µsal
Modelo de colas
Jorge L. Vega Valle
Parámetros del
Intensidad
Tiempo
0.4
3.6412
modelo
2
5
G/G/1
4
5
G/G/s
2
5
4
5
1
=2s
0.4
5.5580
G/G/1/K
= 10 K
0.4
3.6974
G/G/s/K
s = 2,= 10 K
0.4
5.6480
11
G/G/1/∞ /H
=H 5
0.40258
5.5760
4
8
G/G/s/∞ /H
0.41424
7.6892
3
5
2
5
G/G/s/∞ /H con Y
repuestos
s = 4,=H5 s = 6, H = 5,=Y 2
0.41892 0
G/G/∞
7.3626 5.5338
Como conclusión, se puede observar que el tipo del modelo de colas tiene una gran influencia en el tiempo de ejecución del algoritmo de simulación. Así se puede observar como el tiempo de ejecución para el modelo G/G/s/∞ /H es más de dos veces superior al tiempo de ejecución del modelo G/G/1. Esto es debido a que las operaciones que se tienen que realizar no son las mismas en todos los modelos (véase en el apartado 2.3.2 las diferencias entre los distintos algoritmos).
117
Capítulo 4: Conclusiones y trabajo futuro “Las cosas, una vez principiadas, ni se han de olvidar ni dejar, hasta ser acabadas, que es nota de poca prudencia muchos actos comenzados y acabado ninguno” (Mateo Alemán, escritor español, siglos XVI y XVII)
Capítulo 4: Conclusiones
Jorge L. Vega Valle
El desarrollo de la aplicación tiene, entre otras, las siguientes consecuencias directas a la hora de resolver un determinado modelo de colas: Disminuir el tiempo dedicado a la resolución del modelo. Evitar el tener que resolver el modelo “a mano”. Evitar la posibilidad de cometer algún tipo de error en los cálculos. Posibilitar la rápida resolución de otros modelos parecidos. Permitir que se pueda analizar la influencia de la variación de algunos parámetros de entrada en los valores de salida. Ayudar a la toma de decisiones en el caso de que se trate de un modelo real. Permitir resolver por simulación muchos modelos que no tienen una solución analítica. Fomentar la enseñanza y el aprendizaje de la teoría de colas en un ámbito docente.
Es importante resaltar que los alumnos de la Facultad de Informática de A Coruña que han cursado durante este año la asignatura de Teoría de Colas, han estado utilizando esta aplicación para la resolución de modelos reales.
Trabajo futuro
A pesar de que la herramienta es capaz de solucionar los objetivos (ver 1.1.2) que se marcaron al inicio del proyecto, existen otras tareas que se podrían realizar en futuros trabajos: Conseguir dar respuesta a un mayor número de modelos de colas (tanto de la forma M/M/- ó G/G/-, como de redes de colas).
119
Capítulo 4: Conclusiones
Jorge L. Vega Valle
Permitir la posibilidad de ejecutar la aplicación sin la necesidad de tener instalado MATLAB® en el equipo. Sobre este aspecto, caben dos posibilidades: La primera y menos interesante, que es la generación de un ejecutable (por ejemplo, directamente con MATLAB® usando el comando mcc). Para más información sobre como generar un archivo ejecutable con MATLAB®, consultar [7]. La segunda y más interesante, sería el desarrollo de una herramienta de este tipo en un entorno tipo Web (por ejemplo), que permita a distintos usuarios, con la sola presencia de un navegador, poder ejecutar la aplicación. Mejorar la interfaz de la aplicación, de forma que sea más amigable con el usuario.
120
Apéndice A: Características de las de distribuciones usadas en la aplicación A continuación se muestran las características más importantes (función de densidad, media y varianza) de las distintas distribuciones que se utilizan en la aplicación. Distribución exponencial: Exp( ) f ( x) = λ ⋅ e − λ ⋅ x
E( X ) =
Var ( X ) =
x>0
1 λ
1 λ2
Distribución uniforme: U(a,b)
f ( x) =
1
x ∈ (a, b)
b−a
E( X ) =
Var ( X ) =
a+b
2 (b − a ) 2 12
121
Distribución determinista: D(d) f ( x) = d
E( X ) = d
Var ( X ) = 0
Distribución gamma:
(a,p)
f ( x) =
ap ⋅ e − a ⋅ x ⋅ x p −1 Γ( p )
E( X ) =
Var ( X ) =
x>0
p a
p a2
Distribución beta: (p,q,k) x k f ( x) =
p −1
x ⋅ 1 − k
q −1
0
k ⋅ β ( p, q)
k⋅p
E( X ) =
p+q Var ( X ) = k 2 ⋅
p⋅q
( p + q) 2 ⋅ ( p + q + 1)
122
Distribución lognormal: L( , )
f ( x) =
1 x ⋅σ ⋅ 2 ⋅ π
⋅e
E( X ) = e
− (ln( x ) − µ ) 2 2 ⋅ σ2
x∈ R
(µ + σ 2 ) 2
2
2
Var ( X ) = e σ ⋅ e 2 ⋅ µ ⋅ (eσ − 1)
Distribución normal: N( , )
f ( x) =
1 ⋅e σ ⋅ 2 ⋅π
−( x − µ )2 2 ⋅ σ2
x∈R
E( X ) = Var ( X ) = σ 2
Distribución de Weibull: W( , ) f ( x) = λα ⋅ α ⋅ x α −1 ⋅ e − ( λ ⋅ x )
E( X ) =
Var ( X ) =
1 λ
⋅ Γ (1 +
1 α
α
x≥0
)
1 2
λ
2 1 ⋅ Γ(1 + ) − Γ 2 (1 + ) α b
∞
Donde
Γ( p) = ∫ x p −1 ⋅ e − x dx
Además, si p ∈ Ν entonces Γ(p) = (p - 1)!
0
123
Apéndice B: Instalación y ejecución de la aplicación Requisitos de la aplicación Requisitos de software
Para poder ejecutar correctamente la aplicación es necesario tener instalado en el equipo los siguientes componentes: -
Versión de MATLAB® posterior a la 5.0 para así poder visualizar la interfaz.
-
La librería específica que incluye MATLAB® para el cálculo estadístico (Statistics Toolbox), ya que se utilizan funciones de dicha librería para generar números según una determinada distribución (ver 2.1.4 y [2]) Requisitos de hardware
-
Espacio en disco libre: La aplicación necesita un espacio en disco duro ligeramente inferior a 1 Mb y el espacio requerido por MATLAB ® depende tanto de la versión como de los componentes adicionales que instalemos (varía entre 100 Mb y 1.5 Gb)
-
Memoria RAM del equipo: Mayor o igual a 128 Mb.
-
Velocidad del procesador: Mayor o igual a 500 MHz.24
24
Tanto el valor de la memoria RAM, como de la velocidad del procesador, no se corresponden realmente con una limitación, sino más bien con una recomendación para poder ejecutar la aplicación de forma rápida.
124
Instalación de la aplicación
Una vez que tenemos instalado MATLAB® en nuestro equipo, la instalación de la aplicación resulta muy sencilla, pues lo único que tenemos que hacer es indicarle a MATLAB® donde tenemos localizados los archivos que forman parte de la aplicación. Para ello hay que seleccionar en el menú File la opción “Set Path”, y nos aparece una ventana como la que se muestra a continuación:
Figura 11- Localización de la carpeta en donde se encuentra la aplicación
Para especificar el lugar en donde se encuentran instalados los archivos de la aplicación seleccionamos la opción “Add Folder”. Una vez seleccionada la ruta,
125
debemos almacenarla (opción Save), para así no tener que introducirla cada vez que volvamos a ejecutar MATLAB®. Ejecución de la aplicación
Si ya le hemos indicado al sistema donde localizar los archivos que forman parte de la aplicación, estamos en condiciones de poder ejecutarla. Para ello, lo único que hay que hacer es teclear lo siguiente en la Ventana de Comandos de MATLAB®: >> aquas
Figura 12- Ventana inicial de la aplicación
126
Apéndice C: Contenido del CD El contenido del CD adjunto a esta memoria se divide en las siguientes carpetas: Carpeta “Memoria”. Incluye la versión en formato electrónico de este documento. Carpeta “AQUAS”. Incluye los archivos que forman parte de la aplicación. ® Carpeta “Ayuda Matlab”. Incluye los documentos relativos a MATLAB que
han sido consultados con mayor frecuencia. A su vez esta carpeta se divide en las siguientes carpetas: Carpeta “Comenzar con Matlab”. Incluye un archivo [4] en donde se explica lo básico para comenzar a trabajar con MATLAB®. Carpeta “Graficas”. Incluye un documento [5] en donde se explica cuales son los diferentes tipos de gráficas que hay en MATLAB®, así como sus opciones. Carpeta “GUIs”. Incluye un archivo [3] en donde se explica qué hay que hacer para desarrollar una aplicación que trabaje con interfaces usando MATLAB® así como las opciones que proporciona. Carpeta “Manual de Referencia”. Incluye tres documentos [6] en donde se muestra una breve pero muy útil descripción de todas las funciones existentes en MATLAB®. Carpeta “Statistics Toolbox”. Incluye un archivo [2] en donde podemos encontrar la definición de todas las funciones de la librería estadística.
127
Bibliografía Los libros son, de entre mis consejeros, los que más me agradan, porque ni el temor ni la ambición les impiden decirme lo que debo hacer. (Alfonso II de Aragón, el Casto, Rey de la corona de Aragón, siglo XII)
[All-90]
Allen, A.O. (1990). Probability, Statistics and Queueing Theory with Computer Science Applications. Academic Press.
[Bos-02]
Bose,
S.K.
(2002).
An
Introduction
to
Queueing
Systems.
Kluwer
Academic/Plenum.
[Cao-98] Cao, R. y otros (1998). Estadística Básica Aplicada. Tórculo.
[Cao-02] Cao, R. (2002). Introducción a la Simulación y a la Teoría de Colas. Netbiblo. [Dev-86] Devroye, L. (1986). Non-uniform random variate generation. Springer-Verlag.
[Gar-98]
García, A. y otros (1998). Cálculo I. Teoría y problemas de Análisis Matemático en una Variable. Clagsa.
[Gro-98]
Gross, D. and Harris, C.M. (1998). Fundamentals of Queueing Theory. Wiley.
[Hil-97]
Hillier, F. y Lieberman, G. (1997). Introducción a la Investigación de Operaciones. McGraw-Hill.
[Law-91] Law, A.M. and Kelton, W.D.(1991). Simulation Modeling and Analysis. McGrawHill.
128
[Mar-99]
Marchand P. (1999). Graphics and GUIs with MATLAB. CRC Press.
[Pre-97]
Pressman R. G. (1997). Ingeniería del Software. Un enfoque práctico. McGrawHill.
[Red-98] Redfern D. and Campbell C. (1998). The MATLAB 5 Handbook. Springer.
[Wat-89] Watson H.J. and Blackstone J.H. (1989). Computer Simulation. Wiley.
[Wei-95]
Weiss M.A. (1995). Estructuras de Datos y Algoritmos. Addison-Wesley Iberoamericana.
Sitios Web [1] MATLAB. The Language of Technical Computing. Home Page. The MathWorks. http://www.mathworks.com, 2004 [2] MATLAB. The Language of Technical Computing. Statistics Toolbox. The MathWorks. http://www.mathworks.com/access/helpdesk/help/pdf_doc/stats/stats.pdf, 2004 [3] MATLAB. The Language of Technical Computing. Creating Grafical User Interfaces. Version 6. The MathWorks http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/buildgui.pdf, 2004 [4] MATLAB. The Language of Technical Computing. Getting Started With MATLAB . Version 6. The MathWorks http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/getstart.pdf, 2004
129
[5] MATLAB. The Language of Technical Computing. Using MATLAB Graphics. Version 6. The MathWorks http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/graphg.pdf, 2004 [6] MATLAB. The Language of Technical Computing. MATLAB Function Reference. Version 6. The MathWorks http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/refbook.pdf, http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/refbook2.pdf, http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/refbook3.pdf, 2004 [7] MATLAB. The Language of Technical Computing. MATLAB Compiler. Version 6. The MathWorks http://www.mathworks.com/access/helpdesk/help/pdf_doc/compiler/compiler3.pdf, 2004
[8] Manual de Matlab 6.0: “Aprenda Matlab 6.0 como si estuviera en primero”. Escuela Superior de Ingenieros Industriales. Universidad de Navarra, http://www1.ceit.es/asignaturas/Informat1/AyudaInf/aprendainf/matlab60/matlab60.pdf, Julio 2001
130