Diseño óptimo de redes de distribución de agua mediante algoritmos de búsqueda tabú

i Trabajo Fin de Máster Máster Universitario en Organización Industrial y Gestión de Empresas Diseño óptimo de redes de distribución de agua mediante algoritmos de búsqueda tabú Autor: Alicia Robles Velasco Tutor: Pablo Cortés Achedad Dep. Organización Industrial y Gestión de Empresas II Escuela Técnica Superior de Ingeniería Sevilla, 2018 ii iii Trabajo Fin de Máster Máster Universitario en Organización Industrial y Gestión de Empresas Diseño óptimo de redes de distribución de agua mediante algoritmos de búsqueda tabú Autor: Alicia Robles Velasco Tutor: Pablo Cortés Achedad Dep. Organización Industrial y Gestión de Empresas II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 iv v Trabajo Fin de Máster: Diseño óptimo de redes de distribución de agua mediante algoritmos de búsqueda tabú Autor: Alicia Robles Velasco Tutor: Pablo Cortés Achedad El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2018 El Secretario del Tribunal vi vii Agradecimientos A mi tutor, Pablo, y a todos los que enseñando transmiten el gusto por aprender. A mis compañeros de desayuno, porque gracias a ellos he aprendido y crecido cada día. A mi familia, por su apoyo incondicional. Y en particular, a Rafa, por ofrecerme su tiempo y dedicación siempre. “El conocimiento es emocionante, compartirlo y divulgarlo es maravilloso” Francisco Robles Maldonado viii ix Resumen En este Trabajo Fin de Máster se desarrolla un método de optimización que permite resolver el problema de diseño y planificación de redes hidráulicas de distribución. Las redes de abastecimiento son infraestructuras de gran importancia que suponen grandes inversiones. Con el crecimiento demográfico de las ciudades, y el aumento de la contaminación, están apareciendo nuevos retos de diseño y gestión de las mismas. Debido a las características intrínsecas de los sistemas hidráulicos y a la naturaleza discreta de algunas de sus variables, éste es un problema no-lineal mixto y entero. Ante este tipo de problemas de alta complejidad han surgido técnicas de resolución, como las metaheurísticas, que permiten obtener resultados satisfactorios en tiempos reducidos. En este trabajo se diseña un algoritmo de búsqueda tabú ajustando los parámetros a la problemática de estudio. Para la validación del algoritmo se ha resuelto el problema de diseño de redes de distribución aplicado a la red de Alperovits. Esta red ha sido utilizada por muchos autores en la literatura, permitiendo la comparación de los resultados obtenidos en este trabajo con aquellos que obtuvieron otros autores aplicando otras técnicas de resolución. La ejecución del algoritmo se realiza usando el lenguaje de programación Python y, en concreto, una de sus librerías, WNTR, que permite la generación y simulación de redes hidráulicas. xvi Índice de figuras Figura 1 - Distribución del agua del planeta .............................................................................. 3 Figura 2 - Consumo de agua en España y Andalucía ................................................................. 4 Figura 3 - Canal escavado en la roca - Ciudad de Termes, Soria .............................................. 5 Figura 4 - Red de distribución Dos Hermanas, Sevilla, 2018..................................................... 5 Figura 5 - Red ramificada ............................................................................................................. 8 Figura 6 - Red mallada ................................................................................................................. 9 Figura 7 - Red mixta ..................................................................................................................... 9 Figura 8 - Tuberías de fundición dúctil ..................................................................................... 10 Figura 9 - Uniones ...................................................................................................................... 11 Figura 10 - Válvula reguladora .................................................................................................. 11 Figura 11 - Depósitos elevados ................................................................................................. 12 Figura 12 - Disñeo de la Red de Alperovits .............................................................................. 16 Figura 13 - Diseño de Red de Hanoi ......................................................................................... 18 Figura 14 - Diseño de red de New York .................................................................................... 20 Figura 15 - Diagrama de Moody ............................................................................................... 27 Figura 16 - Variables secundarias– Red de Alperovits ............................................................ 31 Figura 17 - Representación de búsqueda del mínimo global ................................................. 35 Figura 18 - Entradas y Salidas de la Búsqueda Tabú ............................................................... 36 Figura 19 - Diagrama de flujo - Búsqueda Tabú....................................................................... 37 Figura 20 - Mecanismo de generación de vecindario ............................................................. 39 Figura 21 - Mecanismo generación de nuevas soluciones ..................................................... 40 Figura 22 - Logos software Python - IDE Pycharm ................................................................... 41 Figura 23 - Estructura del paquete WNTR de Python ............................................................. 37 Figura 24 - Definición de funciones en Python ........................................................................ 40 Figura 25 - Función Alperovits_Network ................................................................................. 40 Figura 26 - Función Decrease_diameter .................................................................................. 41 xvii Figura 27 - Función Change_Network_diameter .................................................................... 41 Figura 28 - Comando para asignar nuevo diámetro ................................................................ 42 Figura 29 - Función Eval_Network_Technical .......................................................................... 42 Figura 30 - Función Invalid ........................................................................................................ 42 Figura 31 - Función Eval_Network_Cost .................................................................................. 43 Figura 32 - Función Tabu_Search .............................................................................................. 43 Figura 33 - Gráficas generadas con la función Graph_solution .............................................. 44 Figura 34 - Gráficas generadas con la función Graph_network ............................................. 45 Figura 35 - Función principal – Python ..................................................................................... 45 Figura 36 - Coste mejor solución vs Número de iteraciones para Tabu-size=4 .................... 58 Figura 37 - Coste mejor solución vs Número de iteraciones para Tabu-size=2 .................... 58 Figura 38 - Evolución Solución actual – Simulaciones solución inicial A ................................ 60 Figura 39 - Evolución Mejor solución – Simulaciones solución inicial A ................................ 61 Figura 40 - Evolución Solución actual – Simulaciones solución inicial B ................................ 63 Figura 41 - Evolución Mejor solución – Simulaciones solución inicial B ................................ 64 Figura 42 - Demanda de los nodos - Simulación 1 .................................................................. 66 Figura 43 - Presión en los nodos - Simulación 1 ...................................................................... 66 Figura 44 - Caudal en los arcos de la red - Simulación 1 ......................................................... 67 Figura 45 - Diámetro de los arcos de la red ............................................................................. 68 Figura 46 - Caudal de los arcos de la red .................................................................................. 69 Figura 47 - Velocidad en los arcos de la red ............................................................................. 70 Figura 48 - Presión en los nodos de la red ............................................................................... 70 Figura 49 - Parámetros de red para la mejor solución encontrada ....................................... 73 1 1. Objeto del Trabajo El presente trabajo fin de máster tiene como objetivo la aplicación de algoritmos de búsqueda tabú al problema de diseño de redes hidráulicas de distribución y la comparación de los resultados con los obtenidos por otros autores usando otras técnicas de optimización. El problema de diseño y planificación de redes hidráulicas de distribución ha sido abordado por diferentes autores en la literatura. No obstante, su estudio no ha sido tan amplio como el de otros problemas de redes, por ejemplo, los de redes de transporte o redes eléctricas. Por ello, parece interesante profundizar en la resolución de este tipo de redes, aportando diferentes estrategias de búsqueda de soluciones. Debido a sus características, el problema de diseño de redes hidráulicas es un problema no lineal mixto y entero. La no linealidad se debe a las ecuaciones que rigen el comportamiento hidráulico del sistema. Además, los diámetros son variables discretas que solo pueden tomar ciertos valores pertenecientes a una gama de diámetros comerciales. Estas características matemáticas hacen que su resolución sea muy compleja, siendo las metaheurísticas la mejor opción para resolver este tipo de problemas encontrada hasta la actualidad. Las metaheurísticas son técnicas capaces de resolver problemas de gran complejidad alcanzando soluciones aproximadas en tiempo reducido. Estos métodos imitan algún comportamiento de la naturaleza para explorar el espacio de soluciones, en general, a través de procedimientos iterativos. La escogida en este trabajo es la búsqueda tabú, la cual simula el comportamiento de la mente humana explorando las posibles soluciones del problema. Esta metaheurística ha sido utilizada para la resolución de otros problemas dando muy buenos resultados, por ello se espera un compartimiento similar para la problemática estudiada en este trabajo. El software utilizado para resolver el problema es Python, un software libre de gran capacidad. Codificar el problema y su correspondiente algoritmo de resolución en lenguaje Python se plantea como una aportación valiosa a la literatura científica relacionada con las redes de distribución de agua. Para el desarrollo de este trabajo fin de máster, en la sección dos se realiza una introducción al tema de estudio. A continuación, en la sección tres, se describe el problema, explicando los tipos de redes de distribución que existen, sus componentes principales y su funcionamiento. En el capítulo cuatro se hace un recorrido por la literatura científica que ha abordado este problema anteriormente para así conseguir una visión del estado actual del mismo, y de los aspectos que no han sido tratados con tanto detalle. Para ello, se analizan las redes modelos más usadas por los investigadores, los modelos matemáticos y sus características, y las técnicas de resolución aplicadas a esta 2 problemática. En el quinto capítulo se presenta la particularización del caso de estudio resuelto en este trabajo, la red modelo simulada, el modelo matemático escogido, el algoritmo de resolución aplicado, y el lenguaje de codificación y software de resolución utilizado. En el capítulo seis se muestran los resultados obtenidos para las diferentes simulaciones realizadas con el fin de validar el modelo. Para ello, se comparan los resultados con los obtenidos por otros autores que resolvieron la misma red hidráulica con distintos métodos. En el capítulo siete se presentan las conclusiones obtenidas a partir de los resultados. Por último, el octavo capítulo recoge las referencias bibliográficas consultadas. Al final del trabajo se añade un anexo con el código completo del algoritmo. 3 2. Introducción El agua ocupa un 70% de la superficie total de la tierra, además es parte constituyente de todos los seres vivos. Del total de agua existente en el planeta se estima que solo el 2.5% se considera agua dulce, siendo el otro 97.5% agua salada. Además, la mayoría del agua dulce, un 69% aproximadamente, se encuentra en la Antártida en estado sólido y otro 30.1% son aguas subterráneas, humedales o acuíferos por lo que se consideran inaccesible para la población. Como consecuencia, el agua dulce, superficial y accesible a los seres humanos se reduce a tan solo el 0.007% del agua que forma el planeta y dicha cantidad se ve reducida cada año debido a la contaminación *1+. Garantizar el acceso a agua potable es uno de los desafíos más importantes del mundo en la actualidad. Figura 1 - Distribución del agua del planeta Fuente: GEO4, 2006 4 En la figura 1 se muestra de manera gráfica la distribución de agua en el planeta en el año 2006. El uso del agua es de lo más variado y dependiendo de este uso el agua requerirá una calidad menor o mayor. Este hecho es importante a la hora de planificar el ciclo del agua en una red, permitiendo el aprovechamiento de aguas para usos que requieran una calidad menor. Existen muchos factores que afectan al consumo de agua como el tamaño de la ciudad o núcleo urbano, el nivel de vida, el nivel de industrialización, el clima, el turismo, etc. En España el consumo de agua se distribuye en: 80% para agricultura o regadío, 14% para abastecimiento de núcleos urbanos y el otro 6% para uso industrial. El abastecimiento de núcleos urbanos incluye el agua apta para el consumo humano, conocida como agua potable, y el agua utilizada para servicios públicos como extinción de incendios. Cada cual cuenta con su red de distribución independiente. Al ser la calidad exigida menor para servicios públicos en muchas ocasiones se utiliza agua recuperada. En las últimas estadísticas del Instituto Nacional de Estadística (INE, 2014) el consumo medio en España era de 132L hab∙día ⁄, cantidad que engloba todos los usos derivados de la vida cotidiana como el aseo, la alimentación, o el uso de electrodomésticos. Estos datos se han utilizado para la confección del siguiente gráfico. Figura 2 - Consumo de agua en España y Andalucía Fuente: Elaboración propia Como se observa en la figura 2, el consumo per cápita de agua en España ronda de los 110 a los 200 litros por habitante y día. Se observa que ha habido una reducción de dicho consumo desde los años 2004-2005, años con consumos más elevados, hasta el 2014, último año del que se tienen datos fiables. Se muestran los datos de Andalucía y España donde puede apreciarse que el consumo de agua en Andalucía se ha reducido bastante con respecto al resto de España en los últimos 20 años. Las redes de abastecimiento surgen ante la necesidad del transporte del agua desde su fuente hasta la población. El origen de las redes de distribución de agua en España se remonta a la civilización romana, ya que está desarrolló infraestructuras de captación, transporte y distribución. Antiguamente las grandes civilizaciones se concentraban cerca 5 de grandes ríos para aprovechar el fácil acceso al agua que este recurso natural ofrece. El desarrollo y la evolución de los métodos de captación, tratamiento y transporte de agua han hecho posible la aparición de grandes centros urbanos en las condiciones más inhóspitas como es el caso de Las Vegas, ciudad localizada en el medio de un desierto. Figura 3 - Canal escavado en la roca - Ciudad de Termes, Soria Fuente: https://www.romanorumvita.com/ La figura 3 se muestra un canal de agua romano de la antigua ciudad de Termes, Soria. Se observa como la estructura estaba preparada para estar cubierta en su totalidad. Figura 4 - Red de distribución Dos Hermanas, Sevilla, 2018 Fuente: EMASESA Como se observa en la figura 4, obra de acondicionamiento y mejora de las redes de abastecimiento en Sevilla, actualmente las redes de distribución también están, generalmente, enterradas. 6 No se podría imaginar la vida moderna como la conocemos si no fuera por el estado actual de las redes hidráulicas tanto de abastecimiento como de saneamiento. El impacto generado por estos sistemas abarca más que la facilidad de obtener el agua potable para el uso diario. El mantenimiento de los índices de calidad y la abundante oferta de este recurso garantizado por las redes de abastecimiento tienen un profundo impacto en la salud de la población. El acceso a agua limpia se relaciona directamente con la reducción de la mortalidad infantil y la prevención de graves enfermedades, entre otras cosas. Hoy en día la población se concentra en núcleos urbanos cada vez mayores, este crecimiento de las ciudades genera la necesidad de ampliar constantemente las redes de abastecimiento. El reto del diseño y la gestión de dichos sistemas es cada días más difícil. El diseño y la construcción de una red hidráulica es un trabajo complejo y ostentoso. Una planificación ineficiente puede suponer grandes gastos públicos o de las empresas encargadas de la construcción y mantenimiento de la red. Por ello, este campo tiene una gran proyección de futuro, cobrando mayor importancia a medida que aparecen nuevos avances tecnológicos que permiten la mejora de estos sistemas. 7 3. Descripción del problema La red de abastecimiento es la encargada de suministrar agua a los consumidores en condiciones adecuadas dependiendo del uso al que esté destinada. El abastecimiento de agua consta de varias fases: Captación, transporte, tratamiento, almacenamiento y distribución. La captación es la fase inicial y consiste en la toma de agua de las diferentes fuentes que existen en la naturaleza. A continuación, el agua es transportada desde las fuentes de abastecimiento hasta los núcleos urbanos donde es tratada y almacenada para su posterior distribución. El tratamiento del agua tiene como objetivo principal conseguir la calidad adecuada del agua dependiendo del uso al que esté destinada, además, se busca mantener las instalaciones de almacenamiento y distribución posteriores en buen estado. El almacenamiento sirve para asegurar el suministro de agua a los consumidores ante diferentes situaciones y regular caudales y presiones en la red de distribución. La última etapa de una red de abastecimiento es la distribución. Este trabajo se centra en dicha etapa que comienza a la salida de la planta de tratamiento de agua o depósitos de almacenamiento, y llega a su fin en el punto de conexión con la instalación de suministro del consumidor. 3.1. Características de las redes de distribución Una red de distribución se define como el conjunto de tuberías y elementos de control que permite el acceso al agua potable a los consumidores. En general, la red de distribución de una ciudad está constituida por multitud de tuberías enlazadas entre sí. Una clasificación de los distintos tramos de la red según su rango podría ser *5+:  Red de transporte: tramo constituido por las tuberías que conectan la planta de tratamiento, depósitos o estación de bombeo con la red arterial. Se caracteriza por contar con tuberías de diámetros mayores.  Red arterial: constituida por los tramos que enlazan los diferentes sectores del área a abastecer.  Red secundaria: tuberías que conectan la red arterial con las acometidas de suministros, tomas contra incendios y bocas de riego.  Acometida: tramo que conecta las instalaciones de suministro directo al consumidor con la red secundaria. Existen muchas clasificaciones de las redes atendiendo a diferentes criterios. Según el tipo 8 de consumo se pueden diferenciar redes para consumo humano, para servicios públicos como extinción de incendios o riego de parques y zonas comunes, y redes para usos industriales. Según el número de depósitos de la red, éstas pueden ser redes de un solo depósito, sin depósitos o con varios depósitos. Según la presión existen redes de alta presión y redes de baja presión. Otra clasificación es entre redes que cuente o no con sistemas de impulsión. Dependiendo de la topología del terreno existen redes en un piso o redes escalonadas. A continuación, se muestran los tipos de redes según su estructura en planta: Ramificada Son redes que cuentan con una tubería principal o primaria de donde parten otras tuberías secundarias y de éstas otras, etc. Normalmente cuentan con un solo nodo fuente situado al comienzo de la tubería principal. El agua circula por la red en un único sentido y se localizan finales de red. Su diseño es sencillo y la regulación de presiones en la misma también, no obstante, en el caso de surgir ampliaciones de la red pueden aparecer presiones ineficientes. En este tipo de redes hay que asegurar la calidad del agua al final del ramal, pues existe mayor tendencia al estancamiento si se produce un bajo uso de algunos ramales. Figura 5 - Red ramificada Fuente: Navarro & Prieto, 2011 En la figura 5 se muestra el diseño de una red ramificada dónde puede observarse la rama principal alimentada por un solo depósito. Mallada En este tipo de redes cada punto de demanda puede estar alimentado por más de una tubería. Generalmente las redes malladas cuentan con más de un depósito, al menos uno de cabecera, en el comienzo de la red y otro en el fin para mantener el equilibrio de presiones en los nodos. A diferencia de las redes ramificadas, en cada tramo el agua puede circular en ambos sentidos. La pérdida de calidad del agua es más improbable, ya que en general, el agua está en circulación. Sin embargo, el diseño de estas redes es más 15 4. Estado del arte En este apartado se realiza un análisis del estado actual del diseño de redes hidráulicas de distribución. Se comienza explicando las redes modelos más utilizadas en la literatura. A continuación, se resumen los modelos matemáticos que representan el funcionamiento de una red hidráulica. Y por último, los enfoques de los diferentes autores para resolver este problema. 4.1. Redes modelo Las siguientes redes que se muestran son las tres más utilizadas en la literatura para simular el comportamiento de algoritmos de resolución aplicados al problema de diseño de redes de distribución de agua potable. Los distintos autores que se enfrentan a este problema y lo hacen con diferentes métodos comparan los resultados obtenidos con los que previamente obtuvieron otros autores. De esta forma es fácil determinar cuan eficaz es un nuevo método de resolución. 4.1.1. La Red de Alperovits and Shamir La primera red que se expone es la red de Alperovits and Shamir (1977). Ésta es una red mallada constituida por 2 mallas, 7 nodos y 8 arcos. El suministro de agua es realizado en totalidad por un único nodo fuente. Es una red de circulación por gravedad, es decir, la red no requiere impulsión. 16 Figura 12 - Disñeo de la Red de Alperovits Fuente: Alperovits and Shamir, 1997 Como se observa en la figura 12 la red está definida por la altura de sus nodos (𝐸𝑖,𝐻𝑖) y la demanda de los mismos (𝑑𝑖). Al ser una red mallada hay que definir hipotéticamente el sentido del flujo de agua (𝑄𝑖,𝑗) y si finalmente, este tiene signo negativo, el sentido será el opuesto. El nodo fuente, nodo 1, es un embalse con una cota constante de 210 metros sobre el nivel del mar. Tabla 1 - Datos referentes a los nodos - Red de Alperovits Fuente: Alperovits and Shamir, 1997 Nodo Altura [m] Mínima presión permitida [m.c.a] Demanda [m3/h] 1 210,0 0,0 -1120,0 2 150,0 30,0 100,0 3 160,0 30,0 100,0 4 155,0 30,0 120,0 5 150,0 30,0 270,0 6 165,0 30,0 330,0 7 160,0 30,0 200,0 17 La tabla 1 recoge la altura de los nodos consumidores de la red y la mínima presión de servicio que se debe satisfacer. Además, la demanda en cada nodo es conocida. Tabla 2 - Longitud de los arcos - Red de Alperovits Fuente: Alperovits and Shamir, 1997 Arco 1 2 3 4 5 6 7 8 Longitud [m] 1000 1000 1000 1000 1000 1000 1000 1000 En la tabla 2 se observa que la longitud de todos los arcos de la red es igual a 1000 metros. Existe un total de 14 diámetros comerciales, cada uno con un coste unitario por unidad de longitud. En la siguiente tabla se recogen dichos diámetros medidos en pulgadas (in) unidades originales, su equivalencia en metros (m), unidades del sistema internacional, y por último el coste unitario asignado a cada diámetro. 1𝑚=39,3701𝑖𝑛 Tabla 3 - Gama de diámetros comerciales - Red de Alperovits Fuente: Alperovits and Shamir, 1997 Diámetros [in] 1 2 3 4 6 8 10 Diámetros [m] 0,0254 0,0508 0,0762 0,1016 0,1524 0,2032 0,254 Coste unitario (u.m./m.L.) 2 5 8 11 16 23 32 Diámetros [in] 12 14 16 18 20 22 24 Diámetros [m] 0,3048 0,3556 0,4064 0,4572 0,508 0,5588 0,6096 Coste unitario (u.m./m.L.) 50 60 90 130 170 300 550 En la tabla 3 se muestra el valor de cada diámetro existente y su coste unitario por unidad de longitud asociado. El número de posibles combinaciones para una solución que represente el diámetro a instalar en cada tubería es 148=1,475,789,056, sin tener en cuenta requisitos técnicos de la red. 4.1.2. Red de Hanoi La red de Fujiwara y Khang (1990), más conocida como red de Hanoi, representa la red de abastecimiento de la ciudad de Hanoi, Vietnam. Es una red mixta formada por 32 nodos y 34 arcos, la distribución de los mismos se observa en la siguiente figura, observándose tres mallas y dos ramas independientes. Al igual que en la red anterior, existe un único nodo fuente, el nodo 1, y es una red de circulación por gravedad, ya que este nodo fuente está a 18 una cota de 100 metros sobre el nivel del mar. Figura 13 - Diseño de Red de Hanoi Fuente: Fujiwara y Khang, 1990 La figura 13 muestra el diseño de la red de Hanoi donde puede verse la distribución de los nodos consumidores y las tuberías que los conectan. La red está definida por los caudales que atraviesan cada arco, para mantener el formato seguido al mostrar los datos de la red de Alperovits, se han calculado los caudales de demanda en cada nodo y son los que se muestran en la tabla 4. 19 Tabla 4 - Datos referentes a los nodos - Red de Hanoi Nodo Altura [m] Mínima presión permitida [m.c.a.] Demanda [m3/h] Nodo Altura [m] Mínima presión permitida [m.c.a] Demanda [m3/h] 1 100 0 0,0 17 0 30 865,0 2 0 30 890,0 18 0 30 1345,0 3 0 30 850,0 19 0 30 60,0 4 0 30 130,0 20 0 30 1275,0 5 0 30 725,0 21 0 30 930,0 6 0 30 1005,0 22 0 30 485,0 7 0 30 1350,0 23 0 30 1045,0 8 0 30 550,0 24 0 30 820,0 9 0 30 525,0 25 0 30 170,0 10 0 30 525,0 26 0 30 900,0 11 0 30 500,0 27 0 30 370,0 12 0 30 560,0 28 0 30 290,0 13 0 30 940,0 29 0 30 360,0 14 0 30 615,0 30 0 30 360,0 15 0 30 280,0 31 0 30 105,0 16 0 30 310,0 32 0 30 805,0 La altura del nodo fuente es de 100 metros sobre el nivel del mar, estando el resto de nodos a nivel del mar. La presión mínima exigida para cada nodo consumidor es, al igual que en la red anterior, de 30 m.c.a, excepto para el nodo fuente que no requiere ninguna presión específica. Tabla 5 - Datos referentes a los arcos - Red de Hanoi Arco 1 2 3 4 5 6 7 8 9 10 Longitud [m] 100 1350 900 1150 1450 450 850 850 800 950 Arco 11 12 13 14 15 16 17 18 19 20 Longitud [m] 1200 3500 800 500 550 2730 1750 800 400 2200 Arco 21 22 23 24 25 26 27 28 29 30 Longitud [m] 1500 500 2650 1230 1300 850 300 750 1500 2000 Arco 31 32 33 34 Longitud [m] 1600 150 860 950 20 En la tabla 5 se observan la longitud en metros de todos los tramos representados en la figura 13. Por último, el número de diámetros disponibles en esta red es considerablemente menor al caso anterior, con un total de seis diámetros, lo que conlleva a un total de 634 = 2.865118∙1026 posibles configuraciones. Tabla 6 - Gama de diámetros comerciales – Red de Hanoi Diámetros [in] 12 16 20 24 30 40 Diámetros [m] 0,3048 0,4064 0,508 0,6096 0,762 1,016 Coste unitario [u.m/m.L] 45,7 70,4 98,4 129,3 180,7 278,3 En la tabla 6 se muestran las medidas de los seis diámetros disponibles en pulgadas (in) y su conversión a metros (m) para mantener el formato seguido en la presentación de la red anterior. El coste unitario por unidad de longitud se presenta en la última fila de la tabla 6. 4.1.3. Red de New York La red de Quindry et al. (1981), más conocida como red de New York, simula la distribución de tuberías en la ciudad de Nueva York (E.E.U.U.). Figura 14 - Diseño de red de New York Fuente: Quindry et al., 1981 21 En la figura 14 se muestra la estructura en plante de la red. Es una red usada por un varios autores para comparar sus métodos de resolución. No obstante, se observa que en la mayoría de los trabajos, los autores han utilizado distintas gamas comerciales de diámetros, por lo que resulta más difícil la comparación cuantitativa de los resultados. Tabla 7 - Datos referentes a los nodos - Red de New York Nodo Mínima altura [ft] Mínima altura [m] Demanda [ft3/s] Demanda [m3/s] 1 300,0 91,4 -2017,5 -57,1 2 255,0 77,7 92,4 2,6 3 255,0 77,7 92,4 2,6 4 255,0 77,7 88,2 2,5 5 255,0 77,7 88,2 2,5 6 255,0 77,7 88,2 2,5 7 255,0 77,7 88,2 2,5 8 255,0 77,7 88,2 2,5 9 255,0 77,7 170,0 4,8 10 255,0 77,7 1,0 0,0 11 255,0 77,7 170,0 4,8 12 255,0 77,7 117,1 3,3 13 255,0 77,7 117,1 3,3 14 255,0 77,7 92,4 2,6 15 255,0 77,7 92,4 2,6 16 255,0 77,7 170,0 4,8 17 255,0 77,7 57,5 1,6 18 255,0 77,7 117,1 3,3 19 255,0 77,7 117,1 3,3 20 255,0 77,7 170,0 4,8 La red está caracterizada por contar con veinte nodos, la mínima altura y la demanda de cada uno de ellos se recoge en la tabla 7. 22 Tabla 8 - Datos referentes a los arcos - Red de New York Arco Longitud [ft] Longitud [m] Arco Longitud [ft] Longitud [m] 1 11600 3535,68 12 12200 3718,56 2 19800 6035,04 13 24100 7345,68 3 7300 2225,04 14 21100 6431,28 4 8300 2529,84 15 15500 4724,4 5 8600 2621,28 16 26400 8046,72 6 19100 5821,68 17 31200 9509,76 7 9600 2926,08 18 24000 7315,2 8 12500 3810 19 14400 4389,12 9 9600 2926,08 20 38400 11704,32 10 11200 3413,76 21 26400 8046,72 La tabla 8 recoge los datos referentes a los arcos de la red, un total de 21 arcos. La red original utiliza los pies (ft) como unidad de medida de longitud. Para mantener el formato de las dos redes modelo mostradas anteriormente, estas medidas se transforman a metros (m), como puede apreciarse en las tablas 7 y 8. 4.2. Modelo de optimización de redes de distribución El diseño óptimo de redes hidráulicas de distribución es un problema en el cual se pretende satisfacer la demanda de agua, con caudal y presión específica, de una serie de nodos consumidores. A lo largo de la historia han aparecido diferentes modelos que simulan el comportamiento de una red hidráulica de distribución, pero todos ellos incluyen las principales características del problema. El punto de partida es la correcta definición de la topología de la red, tanto su estructura como su demanda. La estructura de la red se define por la localización de los nodos, coordenadas y elevación, y la longitud de los arcos que representan las tuberías. La demanda de los nodos consumidores se estima a través de la predicción del crecimiento de la población y el consumo per cápita, datos en general de fácil acceso. El comportamiento de la red se rige por leyes de conservación de masa y energía, ambas se introducen en el modelo en forma de restricción. El objetivo es la elección de un diámetro específico para cada tramo de la red minimizando el coste total a la vez que se satisfacen una serie de restricciones técnicas *10+. Los diámetros asignados a cada tramo de tubería son, en general, las variables principales del modelo, estas variables son discrecionales ya que deben ser escogidas de entre una serie de diámetros comerciales previamente establecidos. Actualmente los modelos más utilizados son los modelos de programación discreta ya que son simples y concisos favoreciendo el alcance de soluciones en relativamente poco tiempo. Un criterio de clasificación de los diferentes modelos matemáticos que simulan el comportamiento de una red hidráulica de distribución puede ser según su función objetivo. Los modelos de redes, en general, buscan minimizar costes, sin embargo, una 23 estrategia común es añadir términos de penalización en la función objetivo por incumplimiento de ciertos requisitos de la red. Por lo que el objetivo del modelo pasa de ser la minimización de los costes de red a incluir la minimización del incumplimiento de las especificaciones. Según este criterio una clasificación de los modelos de redes de abastecimiento podría ser la siguiente *11+:  Modelo sin costes de penalización.  Modelo con costes de penalización por incumplimiento de caída de presiones.  Modelo con costes de penalización por incumplimiento de caída de presiones y velocidad.  Modelo con costes de penalización por incumplimiento de caída de presiones, velocidad y fiabilidad. Cada modelo añade una penalización diferente a la función objetivo. El uso de un modelo u otro dependerá de las características de red y las prioridades de la misma. El modelo general sin costes de penalización es el más usado por los autores que han abordado el problema. Este modelo no contempla términos de penalización en la función objetivo si se incumplen alguna de los requisitos establecidos. Se procede a desarrollar el modelo explicando los conjuntos, los datos y las variables que lo forman. A continuación, se muestra el modelo completo y se realiza una explicación detallada de la función objetivo y de cada una de las restricciones. Conjuntos 𝑁: Conjunto de nodos (tanto consumidores, como nodo fuente) 𝐸𝑖: Conjunto de arcos que llegan al nodo i. ∀𝑖=1,…,𝑁 𝑆𝑖: Conjunto de arcos que salen del nodo i. ∀𝑖=1,…,𝑁 𝐴: Conjunto total de arcos 𝐴=𝐸𝑖∪𝑆𝑖 ∀𝑖=1…𝑁 𝑀: Conjunto de mallas 𝐾: Conjunto de diámetros de la gama comercial 𝑅: Conjunto de números de reales Datos 𝑑𝑖: Demanda del nodo i. ∀𝑖=1,…,𝑁 𝐶𝑜𝑡𝑎𝑖: Cota del nodo i. ∀𝑖=1,…,𝑁 𝐿𝑗: Longitud del tramo o arco j. ∀𝑗=1,…,𝐴 𝐶[∅𝑘]: Coste unitario de tubería para el diámetro ∅𝑘. ∀𝑘=1,…,𝐾 24 Variables ∅𝑘𝑗: Diámetro k ∈ K, gama de diámetros comerciales, asignado al tramo j. 𝑞𝑗: Caudal en el tramo j. 𝑕𝑖: Presión en el nodo i. 𝐻𝑗: Caída de presión en el arco j. Se define el conjunto total de nodos y arcos que forman la red (𝑁 y 𝐴), dos subconjuntos de arcos por cada nodo (𝐸𝑖 y 𝑆𝑖) que representan los arcos que llegan y salen del nodo i. También el conjunto de mallas (𝑀), entendiéndose por malla cualquier conjunto de arcos que forman un recorrido cerrado o bucle en la red. Por último, es necesario definir el conjunto de diámetros posibles a asignar a cada tubería (𝐾). Como se ha comentado, el punto de partida es definir la topología de la red por sus nodos y sus arcos. La posición de los nodos, coordenadas (x,y), es un dato aunque no aparezca en el modelo. 𝐶𝑜𝑡𝑎𝑖 representa la altura sobre el nivel del mar de cada nodo i, y 𝑑𝑖 su demanda expresada en 𝑣𝑜𝑙𝑢𝑚𝑒𝑛 𝑡𝑖𝑒𝑚𝑝𝑜 ⁄. 𝐿𝑗 es la longitud de cada tramo donde debe instalarse una tubería, además se conocen los nodos que une cada uno de los arcos. Cada diámetro de la gama comercial tiene asociado un coste unitario por unidad de longitud 𝐶(∅𝑘). Las variables principales del modelo son los diámetros de las tuberías en cada arco de la red, ∅𝑘𝑗 donde k representa el diámetro de la gama comercial K a instalar en el arco j. Dichos diámetros son variables discrecionales ya que solo pueden tomar ciertos valores predefinidos. La naturaleza de estas variables es una de las causas de la complejidad de la resolución del modelo, explicada en más detalle en apartados posteriores. El resto de variables, o variables secundarias, son los caudales en los arcos, 𝑞𝑗, y las presiones en los nodos, 𝑕𝑖. Estas variables se calculan a través de un proceso iterativo ya que las ecuaciones que las rigen no son lineales. Existen modelos que limitan los valores de velocidad en las tuberías, apareciendo 𝑣𝑗 como variable asociada a cada arco de la red, sin embargo, el modelo tradicional no incluye esta restricción. 31 Figura 16 - Variables secundarias– Red de Alperovits Fuente: Elaboración propia En la figura 16 se han representado los nodos y los arcos que forman la red junto con las variables secundarias del modelo, caudales de agua por los arcos y presiones en los nodos. La tabla 9 recoge todos los datos propios de la red que sirven para caracterizar su estructura y las ecuaciones del modelo, además las unidades de medida de cada dato. La tabla 10 recoge la longitud de los arcos en metros. Tabla 9 - Datos y notación de los nodos - Red de Alperovits Nodo Coordenadas Valor [m] Elevación Valor [m] Demanda Valor [ ⁄] Presión mínima Valor [m.c.a.] 1 ( 1, 1) (2000,2000) 𝐸1 210 𝑑1 -0,311 𝑕1 𝑚𝑖𝑛 - 2 ( 2, 2) (1000,2000) 𝐸2 150 𝑑2 0,028 𝑕2 𝑚𝑖𝑛 30 3 ( 3, 3) (0,2000) 𝐸3 160 𝑑3 0,028 𝑕3 𝑚𝑖𝑛 30 4 ( 4, 4) (1000,1000) 𝐸4 165 𝑑4 0,033 𝑕4 𝑚𝑖𝑛 30 5 ( 5, 5) (0,1000) 𝐸5 155 𝑑5 0,075 𝑕5 𝑚𝑖𝑛 30 6 ( 6, 6) (1000,0) 𝐸6 165 𝑑6 0,092 𝑕6 𝑚𝑖𝑛 30 7 ( 7, 7) (0,0) 𝐸7 160 𝑑7 0,056 𝑕7 𝑚𝑖𝑛 30 32 Tabla 10 - Datos y notación de los arcos - Red de Alperovits Arco Longitud Valor [m] 1 𝐿1 1000 2 𝐿2 1000 3 𝐿3 1000 4 𝐿4 1000 5 𝐿5 1000 6 𝐿6 1000 7 𝐿7 1000 8 𝐿8 1000 Alperovits and Shamir (1997) en su trabajo exponen una gama comercial de 14 de diámetro medidos en pulgadas (in) y su coste unitario por unidad de longitud correspondiente. Se muestra la conversión a unidades del sistema internacional (m) ya que el simulador utilizado para resolver el modelo requiere dichas unidades. Tabla 11 - Datos diámetros comerciales - Red de Alperovits Diámetro Valor [m] (∅ ) [ .𝑚. . .] Diámetro Valor [m] (∅ ) [ .𝑚. . .] ∅1 0,0254 2 ∅8 0,3048 50 ∅2 0,0508 5 ∅9 0,3556 60 ∅3 0,0762 8 ∅1 0,4064 90 ∅4 0,1016 11 ∅11 0,4572 130 ∅5 0,1524 16 ∅12 0,5080 170 ∅6 0,2032 23 ∅13 0,5588 300 ∅7 0,2540 32 ∅14 0,6096 550 En la tabla 11 se recogen el valor y los costes unitarios asociados a cada diámetro comercial. La gama de diámetros obliga a que ningún tramo pueda ser eliminado, como mínimo el diámetro de un tramo será de 0,0254m o lo que es lo mismo 2,54cm. Antes de formular el modelo hay que definir los valores del coeficiente, 𝛼 y la rugosidad 𝐶 usados para el cálculo de las caídas de presiones mediante la ecuación de Hazen-Williams. Estos valores se recogen en la tabla 12. Tabla 12 - Coeficientes de Hazen-Williams para la red de Alperovits Coeficiente Valor 10,667 130 33 El valor del coeficiente es definido en el paquete de simulación hidráulica de Python debido a las unidades que éste utiliza. El coeficiente de rugosidad de Hazen-Williams se fija en el trabajo de Alperovits en 130 *7+. Ambos coeficientes son adimensionales. Modelo 𝐹𝑢𝑛𝑐𝑖ó𝑛 𝑜𝑏𝑗𝑒𝑡𝑖𝑣𝑜: 𝑀𝑖𝑛 [1000∙𝐶(∅ 𝑘1)+𝐶(∅ 𝑘2)+𝐶(∅ 𝑘3)+𝐶(∅ 𝑘4)+𝐶(∅ 𝑘5)+𝐶(∅ 𝑘6) +𝐶(∅ 𝑘7)+𝐶(∅ 𝑘8)] (𝑒𝑐.14) 𝑠𝑎: - 𝐸𝑐𝑢𝑎𝑐𝑖𝑜𝑛𝑒𝑠 𝑑𝑒 𝑐𝑜𝑛𝑠𝑒𝑟𝑣𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑚𝑎𝑠𝑎: 𝑞1=−0,311 (𝑒𝑐.15) 𝑞1−𝑞2−𝑞3=0,028 (𝑒𝑐.16) 𝑞2−𝑞7=0,028 (ec.17) 𝑞3−𝑞4−𝑞5=0,033 (ec.18) 𝑞7+𝑞4−𝑞8=0,075 (ec.19) 𝑞5−𝑞6=0,092 (ec.20) 𝑞6−𝑞8=0,056 (ec.21) - 𝐸𝑐𝑢𝑎𝑐𝑖𝑜𝑛𝑒𝑠 𝑑𝑒 𝑐𝑜𝑛𝑠𝑒𝑟𝑣𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑒𝑛𝑒𝑟𝑔í𝑎: 10,667∙1000 1301,852 ∙*𝑞31,852 ∅𝑘34,87+𝑞41,852 ∅𝑘44,87−𝑞71,852 ∅𝑘74,87−𝑞21,852 ∅𝑘24,87+ =0 (𝑒𝑐.22) 10,667∙1000 1301,852 ∙*𝑞51,852 ∅𝑘54,87+𝑞61,852 ∅𝑘64,87+𝑞81,852 ∅𝑘84,87−𝑞41,852 ∅𝑘44,87+ =0 (𝑒𝑐.23) - 𝑃𝑟𝑒𝑠𝑖𝑜𝑛𝑒𝑠 𝑒𝑛 𝑙𝑜𝑠 𝑛𝑜𝑑𝑜𝑠: 𝑕1=210 (𝑒𝑐.24) 𝑕2=210−10,667∙1000∙𝑞11,852 1301,852∙∅14,87 (𝑒𝑐.25) 𝑕3=𝑕2−10,667∙1000∙𝑞21,852 1301,852∙∅24,87 (𝑒𝑐.26) 𝑕4=𝑕2−10,667∙1000∙𝑞31,852 1301,852∙∅34,87 (𝑒𝑐.27) 34 𝑕5=𝑕4−10,667∙1000∙𝑞41,852 1301,852∙∅44,87 (𝑒𝑐.28) 𝑕6=𝑕4−10,667∙1000∙𝑞51,852 1301,852∙∅54,87 (𝑒𝑐.29) 𝑕7=𝑕6−10,667∙1000∙𝑞61,852 1301,852∙∅64,87 (𝑒𝑐.30) - 𝑅𝑒𝑞𝑢𝑖𝑠𝑖𝑡𝑜𝑠 𝑑𝑒 𝑝𝑟𝑒𝑠𝑖𝑜𝑛𝑒𝑠 𝑒𝑛 𝑙𝑜𝑠 𝑛𝑜𝑑𝑜𝑠: 𝑕2≥30 (𝑒𝑐.31) 𝑕2≥30 (𝑒𝑐.32) 𝑕2≥30 (𝑒𝑐.33) 𝑕2≥30 (𝑒𝑐.34) 𝑕2≥30 (𝑒𝑐.35) 𝑕2≥30 (𝑒𝑐.36) - 𝑉𝑎𝑟𝑖𝑎𝑏𝑙𝑒𝑠: ∅ 𝑘1,∅ 𝑘2,∅ 𝑘3,∅ 𝑘4,∅ 𝑘5,∅ 𝑘6,∅ 𝑘7,∅ 𝑘8 ∈𝐾 𝑞2,𝑞3,𝑞4,𝑞5,𝑞6,𝑞7,𝑞8 ∈𝑅 𝑕2,𝑕3,𝑕4,𝑕5,𝑕6,𝑕7∈𝑅 Hay un total de 21 variables, los 8 diámetros que son variables discretas que deben toman valores de entre un conjunto fijo. Las variables continuas referentes a los caudales en los arcos, en total 7 variables, pues el caudal del arco 1 es conocido e igual al sumatorio de la demanda de los nodos consumidores. Y 6 variables que representan las presiones en los nodos, siendo también variables continuas. La presión en el nodo fuente, 𝑕1, es conocida. La función objetivo depende los diámetros asignados a las distintas secciones o tuberías. Como se muestra en la ecuación 1, esta función es el coste total de la red definido como el sumatorio del coste unitario del tipo de diámetro k asignado al tramo j, 𝐶(∅𝑘𝑗) multiplicado por la longitud de dicho tramo, 𝐿𝑗. Dado, que todos los tramos tienen una longitud de 1000m, este valor sale como factor común en la ecuación. La conservación de masa se impone en las ecuaciones 15-21, existiendo una restricción para cada nodo de la red. Al existir un solo nodo fuente, el caudal del mismo será la suma de las demandas de todos los nodos consumidores, ecuación 15. Para la conservación de energía existe una restricción para cada malla, ecuaciones 22 y 23. Por último, se impone que la presión en cada nodo sea igual o mayor a la mínima presión de servicio establecida, 30 m.c.a., por ello es necesario calcular la presión en cada nodo como la presión en el nodo anterior menos la caída de presión del arco que une dichos nodos. En las ecuaciones 24-30 se calculan las presiones en cada nodo a partir del nodo anterior, no obstante, al ser una red mallada, cada nodo está conectado a más de un arco, 35 por ello, no existe una única forma de calcular dicha presión. En las ecuaciones 24-30 se ha desarrollado el cálculo para uno de los posibles caminos. 5.2. Búsqueda Tabú La Búsqueda Tabú, Tabu Search (TS) en inglés, es una metaheurística que surgió en 1986 por el científico Fred Glover debido a la necesidad de procedimientos rápidos e iterativos para la resolución de problemas complejos de programación entera. El mayor reto al que se enfrentan las metaheurísticas son los óptimos locales, por lo que generar un método capaz de escapar de los mismos representó un gran avance. El algoritmo se basa en el comportamiento de la memoria del ser humano, y especialmente, en su tendencia a desviarse de un curso trazado. Esta tendencia puede verse como un error, pero en muchos casos conlleva a una mejoría. De hecho, es un mecanismo útil, incluso fundamental para el ingenio humano. Partiendo de una solución inicial se construye un vecindario de soluciones adyacentes que pueden ser alcanzadas desde la solución actual. La característica principal de esta metaheurística es la prohibición o penalización de ciertos movimientos, lo que reduce el vecindario en cada iteración. Dichos movimientos son marcados como movimientos “tabú” y permanecerán durante un número de iteraciones con dicho atributo. Una vez ese número de iteraciones, inicialmente fijado, ha transcurrido, el movimiento se desprende del atributo “tabú” y puede ser ejecutado de nuevo sin penalización alguna. Este mecanismo recibe el nombre de memoria a corto plazo y su objetivo es evitar los ciclos en la búsqueda de soluciones. *Glover, 1986+ Figura 17 - Representación de búsqueda del mínimo global Fuente: A. Escudero y P. Aparicio, 2018 En la figura 17 se observa el recorrido que debe seguir la búsqueda tabú para hallar la mejor solución, u óptimo global, del problema. La búsqueda tabú explora el espacio de búsqueda con mecanismos como la memoria flexible, restricciones estratégicas o niveles de aspiración con el fin de escapar de los óptimos locales y alcanzar el óptimo global. Para ello se permite el empeoramiento de la solución como mecanismo para alcanzar una solución final mejor, en la figura 17 se muestra gráficamente este mecanismo. 36 Se definen las principales características del método:  Entorno reducido: al prohibirse ciertos movimientos etiquetados como movimientos tabú, el entorno o vecindario de cada solución se ve reducido y por lo tanto el tiempo total de resolución disminuirá. Es importante fijar el número de iteraciones que un movimiento permanece con el atributo tabú ya que esto influye bastante en los resultados.  Nivel de aspiración: condición que al satisfacerse permite alcanzar una solución, aunque el movimiento que lleva a dicha solución tenga el atributo tabú.  Memoria a corto plazo: almacena características de las soluciones recientes, por ejemplo, el movimiento en la solución anterior que llevó a la solución actual, para asignarle el atributo tabú y así explorar a fondo una sección del espacio de soluciones.  Memoria a largo plazo: el objetivo es intensificar y diversificar la búsqueda almacenando frecuencia de ocurrencias de ciertas características de las soluciones y aplicando criterios para acotar estas ocurrencias. Estas características son las más diferenciadoras de la búsqueda tabú con respecto a otras metaheurísticas. No obstante, el algoritmo puede implementarse incluyendo todas o bien solo alguna de las características anteriores, cada problema tendrá unos requisitos diferentes. A continuación, se procede a explicar mediante un diagrama de flujo el algoritmo de búsqueda tabú diseñado para la resolución de la problemática de estudio. Figura 18 - Entradas y Salidas de la Búsqueda Tabú Fuente: Elaboración propia En primer lugar, el algoritmo requiere dos entradas: la solución inicial y la red de estudio definida, como muestra la figura 18. La solución inicial queda definida por el conjunto de diámetros asignados a cada arco de la red. La red de estudio, en este caso red de Alperovits, se define mediante un paquete de Python, WNRT, explicado con más detalle en la siguiente sección. Una vez aplicada la búsqueda tabú se obtiene la mejor solución, conjunto de diámetros que minimizan los costes de la red, y el coste de dicha mejor solución. A través de este coste se comparan el resto de métodos con el aquí mostrado. 37 Figura 19 - Diagrama de flujo - Búsqueda Tabú Fuente: Elaboración propia 38 En el diagrama de flujo de la figura 19 se muestran los pasos que sigue el algoritmo de búsqueda tabú aplicado al problema en cuestión. Se procede a describir con más detalles dichos pasos:  Definir parámetros - n_iteraciones: es el número de veces que se genera un vecindario. - Tabu_size: es el número de iteraciones que un movimiento permanece en la lista tabú, es decir, con el atributo tabú explicado anteriormente. Ambos parámetros tomarán distintos valores en las diferentes simulaciones con objeto de encontrar una combinación mediante la cual se alcance una solución mejor en menos tiempo.  Inicializar - Sol_actual: la solución actual se inicializa como la solución inicial de entrada al algoritmo. Este parámetro guarda en todo momento la solución actual del algoritmo. Además, en cada iteración se guarda el coste de dicha solución actual en una lista para, una vez realizadas todas las iteraciones, representar la evolución que experimenta. - Mejor_sol: la mejor solución se inicializa como la solución inicial y en cada iteración se actualizará siempre y cuando la solución actual tenga un coste menor que ésta. Por lo que el parámetro Mejor_sol guarda la mejor solución encontrada hasta el momento. Al igual que para la solución actual, el coste de la mejor solución también se almacena en una lista para representar la evolución que experimienta. - Tabu_list: tiene la misma forma que la solución, es decir, si la solución es una lista de 4 elementos la lista tabú será una lista de 4 elementos. Inicialmente, todos los valores de la lista tabú toman el valor de 0, y en cada iteración se incrementará en Tabu_size el elemento que ha sido modificado de la solución en dicha iteración. La lista tabú es el mecanismo de memoria a corto plazo, ya que guarda durante cierto número de iteraciones información de los movimientos anteriores para evitar repetirlos y así no caer en bucles.  Genera_Vecindario La estructura del vecindario es uno de los aspectos más importantes del algoritmo. La generación del vecindario se lleva a cabo a través una función llamada “Reduce_Diámetro” la cual, teniendo como entrada la solución actual, devuelve 39 todas las soluciones originadas al reducir al siguiente diámetro menor de la gama comercial cada uno de los diámetros que forman dicha solución. Figura 20 - Mecanismo de generación de vecindario Fuente: Elaboración propia En la figura 20 se representa a través del ejemplo de una red con cuatro arcos o tuberías, el mecanismo de generación de vecindario. En este ejemplo la gama de diámetros comerciales cuenta con cinco diámetros de dimensiones: 1, 2, 4, 5 y 7 u.L. Se observa que para una solución inicial el vecindario son cuatro soluciones vecinas, en cada una de ellas se ha disminuido un diámetro. Para el caso concreto de la red de Alperovits, al hcontar ésta con ocho arcos, una solución está constituida por una lista de ocho diámetros. El vecindario de dicha solución será un conjunto de ocho soluciones diferentes en las que en cada una se ha reducido uno de los ocho diámetros, y el resto se ha mantenido constante.  Evaluación hidráulica Cada solución vecina se analiza hidráulicamente a través de un simulador del paquete WNTR que reproduce las funciones del software EPANET. Este simulador calcula los caudales en los arcos y las presiones en los nodos aplicando el método de Newton Rapson al sistema de ecuaciones de conservación de masa y energía que rigen el funcionamiento de la red. El requisito para que una solución apruebe la validez hidráulica es que las presiones en todos los nodos, excepto en el nodo fuente, sean mayores a 30 m.c.a.  Genera_Nueva_Solución Este mecanismo sólo se lleva a cabo en el caso de que ninguno de los vecinos de una solución sea hidráulicamente válido. Para escapar de un bucle en el que el algoritmo no avanzaría, se decide aumentar un diámetro, aunque esto suponga un aumento del coste de la solución actual. Para la elección del diámetro a aumentar 40 se escoge al azar entre el conjunto de diámetros cuyo atributo Tabú sea nulo en dicha iteración. Figura 21 - Mecanismo generación de nuevas soluciones Fuente: Elaboración propia En la figura 21 se presenta el mecanismo de generación de nuevas soluciones, para el supuesto de que el vecindario de una solución actual esté vacío. En el ejemplo de la figura, los únicos diámetros con estado tabú igual a 0 son los diámetros de los arcos i igual a 1 e i igual a 3, se escoge aleatoriamente uno de ellos y se aumenta su diámetro. Una vez generada esta nueva solución se actualiza la solución actual, Sol_actual, y la lista tabú evitando que en las próximas iteraciones se vuelva a disminuir el mismo diámetro que se ha aumentado.  Evaluación Coste Se evalúa el coste de las soluciones vecinas que son válidas hidráulicamente. Este coste se calcula como el coste unitario del diámetro escogido para cada arco por su longitud. Una vez calculado dicho valor se pueden comparar soluciones y así poder escoger la de menor coste total. A continuación, se comprueba si la solución con menor coste del vecindario tiene estado tabú. En el caso de que el estado tabú sea 0, se actualiza la Sol_actual, el coste de la solución actual y la lista tabú, aumentando en Tabu_size el elemento que se modifica para alcanzar dicha solución. Si la solución actual tuviese un coste menor al coste de la mejor solución, ésta también se actualizaría. Si por el contrario, la mejor solución del vecindario tiene estado tabú, se aplica el criterio de aspiración. 42 Figura 28 - Comando para asignar nuevo diámetro Fuente: Elaboración propia En la figura 28 se presentan los comandos utilizados para asignar un nuevo diámetro a una tubería de la red “wn”. Con la sentencia get_link se elige un arco, llamado en este ejemplo pipe1. En la siguiente línea se extrae el parámetro “diámetro” de dicho arco y se le asigna el valor de la lista de diámetros D que le corresponda. Figura 29 - Función Eval_Network_Technical Fuente: Elaboración propia La función Eval_Network_Technical, figura 29, simula la red usando el módulo EpanetSimulator del paquete WNTR. En primer lugar, se realiza la simulación hidráulica de la red verificando que existe una solución definida por los caudales en los arcos y las presiones en los nodos. Una vez se alcanza dicha solución se comprueba si la presión en todos los nodos de la red supera la presión mínima de servicio. Las entradas de la función son la red previamente definida, wn, y la lista de diámetros asignados a cada tubería, D. La salida será True si dicho conjunto de diámetros supone una solución válida y la presión en todos los nodos supera la mínima. Si se incumple alguno de estos requisitos la salida será False. Figura 30 - Función Invalid Fuente: Elaboración propia La figura 30 presenta las entradas y salidas de la función Invalid. Esta función se crea para evitar realizar comprobaciones hidráulicas de soluciones ya estudiadas. Es decir, 43 ahorrar tiempo de resolución. La función comprueba si una configuración determinada de diámetros, lista D, ya ha sido catalogada como solución inválida hidráulicamente. Para ello, cada vez que se encuentra una solución inválida, ésta se añade a la lista Invalid_solutions. Invalid_solutions es, por tanto, una lista de listas. En la función se recorre dicha lista comparando las soluciones inválidas almacenadas con la solución de entrada, lista de diámetros D. La salida será True si la lista de diámetros D se corresponde con alguna de las listas contenidas en Invalid_solution, y False, en caso contrario. Figura 31 - Función Eval_Network_Cost Fuente: Elaboración propia La función Eval_Network_Cost, figura 31, calcula el coste de la red para unos diámetros específicos en cada arco. Esta función tiene una sola entrada, la lista de diámetros D, y una sola salida, el coste, en unidades monetarias, de dicha asignación de diámetros. Para realizar este cálculo la función contiene los costes asociados a cada diámetro de la gama comercial, establecidos por Alperovits and Shamir en su trabajo. Se recorre la lista de diámetros D, comparando cada elemento de la lista con los valores de los diámetros de la gama comercial, si coinciden se añade el coste de dicho diámetro multiplicado por la longitud del arco en el que está instalado. Figura 32 - Función Tabu_Search Fuente: Elaboración propia 44 La función Tabu_Search, figura 32, es la principal del código, conteniendo los mecanismos referentes a la metaheurística presentados en la sección 5.2. Dentro de esta función se llama a las otras funciones explicadas, ya que, en general, son parte de los procedimientos incluidos en la búsqueda tabú. A partir de la solución inicial, lista de diámetros D, la función devuelve la mejor solución encontrada tras aplicar la metaheurística, best_solution, y el coste, calculado con la función Eval_Network_Cost, de dicha solución. Por último, se definen dos funciones encargadas de generar distintas gráficas de los resultados obtenidos. Éstas son Graph_solution y Graph_network. A través de la función Graph_solution se representan dos gráficas que muestran la evolución de la mejor solución y la solución actual durante la ejecución de la búsqueda tabú. Las entradas de esta función son dos listas, una que contiene los costes de las mejores soluciones encontradas durante la implementación del algoritmo, y la otra conteniene todas las soluciones que toman el estatus de solución actual. Figura 33 - Gráficas generadas con la función Graph_solution Fuente: Elaboración propia En la figura 33 se exhiben las dos gráficas que se crean al utilizar la función explicada. A la izquierda la evolución de las mejor solución, y a la derecha la evolución de la solución actual. La función Graph_network simula la red generando tres gráficas: la demanda y los caudales en los arcos, y la presión de los nodos. Esto se realiza utilizando el módulo Graphics del paquete WNTR. 45 Figura 34 - Gráficas generadas con la función Graph_network Fuente: Elaboración propia En la figura 34 se observan dos de las gráficas que se generan al aplicar la función Graph_network. Estas gráficas permiten comprobar el rango en el que se mueven los parámetros principales de la red. Por último, en la siguiente figura se muestra la función main(), función principal del código. Figura 35 - Función principal – Python Fuente: Elaboración propia Las líneas de código que forman la función están recogidas en la figura 35. Una de las características de la función main() es que no tiene entradas ni salidas, al ejecutar el código ésta es la función que se ejecuta por defecto. Se comienza estableciendo la solución de partida, o solución inicial, lista de diámetros asignados a cada tramo de la red y se llama a la función Tabu_Search, cuya entrada es dicha lista. Además, se mide el tiempo de ejecución del algoritmo de búsqueda tabú. 5.3.4. Pseudocódigo “Pseudo” deriva del griego “seudo” y su traducción sería falso, se utiliza para referirse a algo que no es original o que es una imitación. “Código” proviene del latín, representando el documento donde se establecían las leyes. En la actualidad, un código es un conjunto de símbolos con un determinado valor, previamente asignado, que se combinan respetando ciertas reglas y permiten transmitir un mensaje. En programación, código se utiliza para referirse al lenguaje que interpretan los procesadores o softwares. Mientras que los 46 pseudocódigos describen dichos algoritmos para una fácil interpretación humana. Para concluir el capítulo 5 se presenta el pseudocódigo del algoritmo implementado. El código completo se añade en un anexo al final del trabajo. ----- INICIALIZACIÓN ----- Sol_Inicial=Definir_Solucion_Inicial() Wn=Definir_red() Inicializar n_iteraciones, Tabu_size, Tabu_list Inicializar Lista Soluciones NO_validas Inicializar Lista_Mejor_solucion, Lista_Solucion_Actual Sol_actual=Sol_Inicial Mejor_solución=Sol_actual Mientras (n<n_iteraciones): ----- VECINDARIO ----- Vecindario=Reduce_diametro(Sol_actual,wn) Quitar del vecindario vecinos ∈ NO_validas Evaluación Hidráulica(Vecindario) Actualizar lista NO_validas si Vecindario Vacío: i =Indice_aleatorio ∈ Lista Tabu=0 Sol_actual=Incrementa_diametro(Sol_actual,i) sino: CosteVecindario=Eval_Coste(Vecindario) vecino=vecino ∈ Vecindario con menor CosteVecindario Coste_vecino=Eval_Coste(vecino) ----- Criterio Aspiración y Memoria corto plazo----- si vecino∉ Lista_tabu or Coste_vecino<Coste_Mejor_sol: Actualiza Sol_actual Coste_Sol_actual=Coste_vecino Añadir Sol_actual a Lista_Solucion_Actual 47 Actualiza Tabu_list incrementando en Tabu_size+1 si Coste_Sol_actual<Coste_Mejor_sol: Actualiza Mejor_sol Coste_Mejor_sol=Coste_Sol_actual Añadir Mejor_sol a Lista_Mejor_Solucion Tabu_list=Tabu_list-1 n=n+1 ----- Gráficas ----- Genera grafica Lista_Mejor_Solucion Genera grafica Lista_Solucion_Actual Genera graficas parámetros Red FIN El pseudocódigo anterior se divide en inicialización, generación del vecindario, criterio de aspiración, memoria a corto plazo, y gráficas. En la inicialización se define la solución inicial y la red, además, se inicializan los parámetros característicos de la búsqueda tabú. La generación del vecindario se realiza a través de los mecanismos explicados en los apartados anteriores. Una vez generado el vecindario, éste se evalúa y se analiza la mejor solución del mismo, si esta solución tiene estatus tabú se aplica el criterio de aspiración. En el caso de que no existe ninguna solución válida en un vecindario, se aplicaría el procedimiento de aumentar uno de los diámetros de la solución actual. Una vez se ha implementado el algoritmo tantas veces como iteraciones se fijen al inicio, se generan las gráficas. Con este pseudocódigo queda explicado el algoritmo de resolución. 48 6. Resultados En esta sección se muestran las simulaciones realizadas y los resultados obtenidos con el fin de validar el algoritmo. En cada simulación se modifican ciertos parámetros que afectan al método de resolución. En primer lugar se explica la elección de dichos parámetros, a continuación se muestra los resultados obtenidos en las simulaciones realizadas. En los siguientes apartados se evalúa el algoritmo y las características hidráulicas de las soluciones alcanzadas. Por último, se comparan los resultados obtenidos con los de otros autores recogidos en la literatura. 6.1. Ajustes de los parámetros del algoritmo En esta sección se estudian los distintos parámetros que afectan al método de resolución con el fin de ajustar sus valores y así evitar ineficiencias computacionales. Los parámetros que más influyen en la ejecución del algoritmo de búsqueda tabú son: la solución inicial, el criterio de parada, y el tabu_size, o número de iteraciones que un movimiento permanece con el estatus “tabú”.  Solución inicial En primer lugar, se ha decidido ejecutar el algoritmo para dos soluciones iniciales diferentes. Esto permitirá comprobar la validez del algoritmo sea cual sea la solución de partida. Las soluciones se caracterizan por una lista con ocho elementos que representan el diámetro asignado a cada arco de la red, por ejemplo, el diámetro del arco 1 de la red es el primer elemento de la lista, y así sucesivamente. La solución inicial comienza tomando el papel de solución actual y , entonces, el algoritmo generará el vecindario de dicha solución. Tabla 14 - Soluciones de partida Solución inicial A *0.6096,0.6096,0.6096,0.6096,0.6096,0.6096,0.6096,0.6096+ Solución inicial B *0.508,0.4064,0.4064,0.254,0.0508,0.3048,0.3048,0.1016+ La primera fila de la tabla 14 muestra la primera solución inicial utilizada, de aquí en adelante nos referiremos a ella como solución inicial A. Ésta es una solución en la cual todas las tuberías de la red tienen el mayor diámetro de la gama comercial. El coste asociado a esta solución es de 4,400,000 u.m. Esta solución es inadmisible hidráulicamente 49 ya que las presiones en los nodos no estarían dentro de los límites establecidos, además, la velocidad del agua en las tuberías sería excesivamente baja debido a los grandes diámetros de las mismas. No obstante, parece lógico como norma general comenzar por la peor solución existe para así explorar un espacio de búsqueda mayor. La segunda solución de partida, solución inicial B, es una solución válida hidráulicamente y con un coste de 498,000 u.m. Se pueden observar los diámetros, en metros, asignados a cada arco de la red en la segunta fila de la tabla 14. Esta solución ya reduce en un 88.68% el coste de la solución inicial A. No obstante, dado que se conocen soluciones aún mejores para este problema, se sabe que ésta no es una solución óptima y por tanto, plantea un escenario diferente al anterior.  Tabu_size Otro de los parámetros que ha sido modificado para las distintas simulaciones es el tabu_size, número de iteraciones que un movimiento permanece en la lista de movimientos prohibidos o lista tabú. Generalmente este parámetro se fija en torno a la mitad del total de movimientos posibles partiendo de una solución, es decir, si el vecindario de una solución son 10 soluciones vecinas, el tabu_size más común es 5. No hay reglas establecidas, pero un número demasiado grande de este parámetro puede generar vecindarios vacíos, lo que haría que el tiempo de resolución aumentase sin avances en el estudio del espacio de soluciones. En cambio, si el tabu_size es demasiado pequeño la metaheurística no cumpliría su función ya que no prohibiría dichos movimientos en un número significativo de iteraciones. Dado que el vecindario de una solución está formado por ocho soluciones vecinas, se elige un tabu_size igual a 4. Hay que tener en cuenta que por las características del problema, existirán soluciones que no sean válidas hidráulicamente y, por ello, no se añadirán al vecindario, apareciendo en ciertas ocasiones vecindarios de menor tamaño. En consecuencia, se decide implementar el algoritmo también utilizando un tabu_size igual a 2. En resumen, se realizaran simulaciones con estos dos valores para el parámetro tabu_size.  Criterio de parada El criterio de parada influye en el tiempo de resolución del problema, tiempo de cálculo del software. Es conveniente comenzar con un número de iteraciones reducido e ir aumentándolo, para comprobar el correcto funcionamiento del algoritmo. Existen dos criterios generales para fijar este parámetro. El primero, y más simple, es fijar un número máximo de iteraciones. En ocasiones se observa que a partir de un cierto número de iteraciones la mejor solución del problema no se reduce. Por ello, el segundo criterio es crear un mecanismo de control del número de iteraciones en las que la mejor solución no se actualiza, y detener el algoritmo cuando el número máximo fijado se alcance. En este trabajo se escoge el primer criterio, estableciendo un número de iteraciones para las diferentes simulaciones de 100, 1000 y 2000 iteraciones. Cabe destacar el algoritmo se testó para 5000 y hasta 10000 iteraciones, no observándose mejora a partir de las 2000 iteraciones. Por ello, este número se establece como el máximo. Un criterio establecido por Baños et al. (2010) en su trabajo, en el cual se resuelven redes de distintos tamaños, es el siguiente: 50 𝑛𝑖=𝐾∗𝑛𝑎∗lo (𝑛∅) (𝑒𝑐.36) Dónde: 𝑛𝑖: Número de iteraciones 𝐾: Valor constante 𝑛𝑎: Número de arcos de la red 𝑛∅: Número de diámetros de la gama comercial En este trabajo al resolver una sola red se decide utilizar un número fijo de iteraciones. No obstante, este criterio tiene utilidad para futuras pruebas del algoritmo en otras redes modelo. Una vez definido el valor de los parámetros que interesa modificar, se obtiene un total de doce combinaciones, por lo que se realizarán doce simulaciones con las siguientes características. Tabla 15 - Parámetros de las distintas simulaciones Tabu_size 4 2 Número de iteraciones 100 Solución inicial A Solución inicial A Solución inicial B Solución inicial B 1000 Solución inicial A Solución inicial A Solución inicial B Solución inicial B 2000 Solución inicial A Solución inicial A Solución inicial B Solución inicial B En la tabla 15 se observa el valor de los parámetros de cada una de las doce simulaciones realizadas. 6.2. Simulaciones Se procede a mostrar en detalle los datos de entrada y los resultados obtenidos en cada una de las doce simulaciones realizadas. Tabla 16 - Entradas simulación 1 Arco 1 2 3 4 5 6 7 8 Diámetro [m] 0,6096 0,6096 0,6096 0,6096 0,6096 0,6096 0,6096 0,6096 Coste de la solución [u.m.] Tabu_size Número de iteraciones 4.400.000 4 100 51 Tabla 17 - Salidas simulación 1 Coste de la solución [u.m.] Tiempo de resolución [s] Nº soluciones no válidas 420.000 35 336 Arco 1 2 3 4 5 6 7 8 Diámetro [m] 0,5080 0,2540 0,4064 0,0254 0,3556 0,2540 0,2540 0,0254 Caudal [m3/s] 0,311 0,102 0,181 0,000 0,147 0,056 0,075 0,000 Velocidad [m/s] 1,535 2,019 1,395 0,541 1,484 1,100 1,471 -0,367 Nodo 1 2 3 4 5 6 7 Presión[m.c.a.] 0,00 55,96 30,87 46,56 32,48 30,80 30,90 Las entradas de la primera simulación, recogidas en la tabla 16, son la solución inicial, diámetros asignados a cada arco de la red, el coste de dicha solución en unidades monetarias, el valor del tabu_size y el número de iteraciones. Las salidas o resultados de la primera simulación, expuestas en la tabla 17, son la mejor solución encontrada y el coste que conlleva. El tiempo de resolución, parámetro que no aporta mucha información ya que depende del procesador y la memoria del ordenador en el que se realizan las simulaciones. No obstante, ya que todas las simulaciones se realizan con el mismo ordenador, Intel Core i5-5200U con 4GB de memoria, el tiempo de resolución sirve para comparar unas simulaciones con otras. Otra de las salidas almacenadas es el número de soluciones no válidas analizadas, este parámetro sólo almacena cada solución no válida hidráulicamente una vez. Permite valorar el entorno analizado en cada simulación. Si el número de soluciones no válidas analizadas es muy alto, se ha dedicado mucho tiempo de cálculo a explorar un espacio de soluciones sin validez hidráulica, lo que no conviene si el objetivo es reducir el tiempo total de resolución. Los parámetros característicos de la red cuya distribución de diámetros es la solución final se obtienen usando el paquete de Python WNTR. Como se observa en la tabla 17, el algoritmo facilita el caudal en cada arco y la presión en cada nodo de la red. Además, la velocidad media del agua en cada arco se calcula a partir del caudal y del diámetro de cada tubería usando la siguiente fórmula. 𝑣𝑗=𝑄𝑗 𝑆𝑗 (𝑒𝑐.37) Dónde: 𝑣: Velocidad media del agua en el arco j 𝑄𝑗: Caudal que circula a través del arco j 𝑆𝑗: Sección transversal de la tubería en el arco j En este trabajo la sección transversal de las tuberías es circular, por ello se calcula siguiendo la siguiente ecuación. 58 Figura 36 - Coste mejor solución vs Número de iteraciones para Tabu-size=4 Fuente: Elaboración propia Figura 37 - Coste mejor solución vs Número de iteraciones para Tabu-size=2 Fuente: Elaboración propia Las gráficas de las figuras 36 y 37 muestran el coste de la mejor solución encontrada para las doce simulaciones. En azul las simulaciones cuya solución de partida era la solución inicial A, y en rojo las simulaciones cuya solución de partida en la solución inicial B. En la figura 36 se observan los resultados obtenidos para un tabu_size igual a 4, mientras que la figura 37 presenta los resultados para un valor del tabu_size igual a 2. Se observa que tanto para un tabu_size igual a 4 como a 2 se obtienen mejores resultados cuando se parte de la solución inicial A, es decir, la peor solución posible del problema. Esto puede ocurrir debido a que la probabilidad de escapar de un óptimo local en pocas iteraciones es menor cuando se parte de una solución buena. Al incluir el algoritmo una generación aleatoria de soluciones cuando no se encuentran soluciones válidas, parece 380000 390000 400000 410000 420000 430000 440000 450000 460000 470000 480000 490000 100 1000 2000 Coste mejor solución [u.m.] Iteraciones Sol Inicial A Sol Inicial B 380000 390000 400000 410000 420000 430000 440000 450000 460000 470000 480000 490000 100 1000 2000 Coste mejor solución [u.m.] Iteraciones Sol Inicial A Sol Inicial B 59 que esta aleatoriedad funciona mejor si partimos de una solución de mala calidad. La función Increase_diameter aumenta un diámetro escogido al azar cuando en un vecindario no se encuentra ninguna solución que cumpla los requisitos marcados. En ambas gráficas se observa que cuando el número de iteraciones aumenta las soluciones encontradas son, en general, mejores. Como se comentó en la anterior sección, se ejecuta el algoritmo para un número mayor de iteraciones sin observarse mejora en la solución, por ello, se decide parar en 2000 iteraciones. 6.3. Evaluación del algoritmo Una forma útil de verificar si el algoritmo de búsqueda tabú está funcionando correctamente es representando las gráficas de evolución, tanto de la mejor solución, como de la solución actual. Una característica distintiva de esta metaheurística es que permite el empeoramiento de la solución con el fin de escapar de los óptimos locales. Se recuerda que la solución actual se actualiza en cada iteración, mientras que la mejor solución solo se actualiza si la solución actual tiene un coste menor al de la mejor solución encontrada. A continuación, se muestran las gráficas de evolución, tanto de la mejor solución, como de la solución actual para las doce simulaciones realizadas. De nuevo el color azul representa las simulaciones cuya solución de partida es la solución inicial A, y el color rojo las que utiliza la solución inicial B. El eje X representa la iteración y el eje Y el coste de la solución para dicha iteración en unidades monetarias. Las gráficas se han organizado en dos columnas, la primera correspondiente a aquellas simulaciones cuyo tabu_size es igual a 4, mientras en la segunda columna el tabu_size es igual a 2. Cada fila de figuras representa un número de iteraciones, mostrado en la parte derecha. 60 Figura 38 - Evolución Solución actual – Simulaciones solución inicial A Fuente: Elaboración propia 61 Figura 39 - Evolución Mejor solución – Simulaciones solución inicial A Fuente: Elaboración propia Las figuras 38 y 39 presentan la evolución de las soluciones para las seis simulaciones cuya solución de partida es la solución inicial A, esta solución se caracteriza por ser la peor del 62 problema. Se observa cómo tanto la solución actual como las mejores soluciones tienen un coste inicial muy alto 4,400,000 u.m., y por lo tanto, se aprecia un descenso del mismo muy acentuado en las primeras iteraciones. El número total de iteraciones de cada simulación puede observarse en el eje X de las gráficas de solución actual, con título Current_Solution. El funcionamiento del algoritmo queda validado al observar que la solución actual empeora y consigue mejorar después, fenómeno que ve más fácilmente en las dos primeras simulaciones de la figura 38. Con este mecanismo se consigue analizar un espacio de soluciones mayor. En las gráficas de evolución de la mejor solución el eje X representa el número de veces que la mejor solución se actualizó, es decir, que se encontró una solución mejor. Se observa que, independientemente del número de iteraciones, 100, 1000 o 2000, la mejor solución se actualiza alrededor de 50 veces. En las primeras iteraciones el algoritmo encuentra mejores soluciones con facilidad, no obstante, cuando el problema alcanza una solución con un coste relativamente bajo es cuando le cuesta más escapar de los óptimos locales. 63 Figura 40 - Evolución Solución actual – Simulaciones solución inicial B Fuente: Elaboración propia 64 Figura 41 - Evolución Mejor solución – Simulaciones solución inicial B Fuente: Elaboración propia En las figuras 40 y 41 se muestran las gráficas de evolución para las simulaciones cuya solución de partida es la solución inicial B, solución cuyo coste es 498,000 u.m. En todas las 65 simulaciones de este conjunto se observa como el coste de la solución actual, Current_solution, experimenta un ascenso considerable en las primeras iteraciones, figura 40. Al partir de una solución ya optimizada, el algoritmo encuentra dificultades para hallar una mejor solución por lo que comienza aumentando ciertos diámetros, incrementando así el coste total de la solución en las primeras iteraciones. Se observa como gracias a este empeoramiento se consigue explorar nuevos espacios de búsqueda alcanzando finalmente mejores soluciones. A diferencia de la evolución de las mejores soluciones para el conjunto de simulaciones anteriores, en estas seis simulaciones la mejor solución se actualiza entre 3 y 14 veces, como muestra la figura 41. Al partir de una solución mejor, el número de veces que esta puede mejorar se reduce. Se puede observar que, para las dos primeras simulaciones de este conjunto, habiendo un total de 100 iteraciones solo en tres ocasiones se encuentra una solución mejor, además, en dichos casos la soluciones alcanzadas no tiene un coste bajo en comparación con el resto. Se concluye que el algoritmo diseñado obtiene mejores resultados utilizando una solución de partida de baja calidad ya que esto permite explorar un espacio de búsqueda mayor. Para este tipo de soluciones de partida, la respuesta mejora usando un tabu_size igual a 4, es decir, a la mitad del tamaño del vecindario. En este punto, queda comprobado el correcto funcionamiento del algoritmo ya que éste explora un espacio de soluciones amplio obteniendo soluciones satisfactorias en tiempo reducido. En el siguiente apartado se analiza la validez hidráulica de las soluciones obtenidas. 6.4. Evaluación hidráulica de las soluciones Para comprobar las características hidráulicas de la red ante las distintas soluciones se han programado una serie de funciones en Python que permiten la representación en gráficas de los parámetros fundamentales. Las gráficas se generan utilizando el paquete Matplotlib de Python, permitiendo controlar de forma visual y rápida, el cumplimiento general de los requisitos técnicos de la red. La presión en los nodos no puede ser menor a 30 m.c.a., excepto para el nodo fuente, en el que no se impone ningún límite. Aunque estas gráficas se extraen de todas las simulaciones realizadas, solo se muestran a continuación, las relativas a la solución encontrada en la primera simulación. 66 Figura 42 - Demanda de los nodos - Simulación 1 Fuente: Elaboración propia Figura 43 - Presión en los nodos - Simulación 1 Fuente: Elaboración propia 67 Figura 44 - Caudal en los arcos de la red - Simulación 1 Fuente: Elaboración propia Las figuras 42, 43 y 44 muestran los parámetros característicos de la red para la mejor solución encontrada en la simulación 1. En la figura 42 se observa la demanda de cada nodo, esto era un dato de partida, no obstante, a través de esta gráfica se comprueba que dicha demanda es atendida. Recordad que las unidades de medida de la demanda son 3 ⁄, por lo que varía de -0.311 3 ⁄ para el nodo 1, o nodo fuente, hasta 0.092 3 ⁄ para el nodo 6, nodo con mayor demanda. Un nodo con demanda negativa es un nodo que suministra agua a red, mientras que una demanda positiva significa que el nodo es consumidor de agua. En la figura 43 se observan las presiones correspondientes a cada nodo de la red. A simple vista se comprueba que todas superan los 30 m.c.a., excepto la presión del nodo 1 que no tenía límite inferior. Se puede comprobar que los nodos 2 y 4 tienen una presión mayor al resto de nodos. Es natural que cuanto más cerca esté un nodo del nodo fuente (nodo 1), mayores presiones alcanza. Por último, la figura 44 permite observar la cantidad de caudal que atraviesa cada arco de la red. Esta gráfica muestra el recorrido del caudal, siendo útil para comprobar que tuberías tendrán un uso mayor. En la red representada se observa que los arcos 4 y 8, transportan significativamente menos caudal que el resto de arcos, ya que son de color azul intenso. Se procede a representar en diferentes gráficas el diámetro, el caudal y la velocidad de cada arco de la red para las distintas soluciones encontradas. Además, se representan las presiones en los nodos. En las doce simulaciones realizadas se encuentran cinco soluciones diferentes, ya que en muchas de las simulaciones se llega a la misma solución final. 74 7. Conclusiones El objetivo de este trabajo fin de máster era el diseño, desarrollo y validación de un algoritmo de búsqueda tabú para la resolución del problema de diseño de redes de distribución de agua. A lo largo del trabajo se ha explicado la problemática y su situación actual. El crecimiento demográfico de las ciudades incurre en constantes ampliaciones de las redes de distribución. Además, el desarrollo de tecnología avanzada y la aparición de nuevos materiales de construcción, junto con el deterioro que sufren las redes, hace que la reposición de tramos, o incluso redes completas, sea una tarea cada día más común. Ante esta situación, cobra gran importancia la mejora del diseño y la planificación de las redes hidráulicas de distribución, problemática que carece de atención en la literatura científica en comparación a otros problemas de redes, como son los de red de transporte o red eléctrica. Para introducir el tema, en este trabajo se definen los tipos de redes, sus componentes y las funciones principales que éstas desarrollan. Una vez entendida la problemática, se explica el modelo matemático tradicional, que incluye las ecuaciones de conservación de masa y energía, propias de las redes hidráulicas. Las características del modelo hacen de él un problema no lineal mixto y entero. Esto se debe a la naturaleza discreta de los diámetros, ya que éstos deben ser escogidos de entre los disponibles en una gama comercial. Y a la no linealidad de las ecuaciones de conservación de masa y energía. El algoritmo de búsqueda tabú diseñado explora el espacio de soluciones del problema, definido por el número de variables del mismo. Para la aplicación del algoritmo a un caso concreto se escoge la red de Alperovits, red famosa en la literatura ya que muchos autores han publicado los resultados obtenidos al aplicar distintas técnicas de resolución a dicha red. Se utiliza el lenguaje Python para programar todos los procedimientos y cálculos de resolución diseñados. Este lenguaje cuenta con multitud de librerías, las cuales sirven de soporte para realizar funciones específicas. Se realizó una búsqueda encontrando una librería llamada WNTR que permite la definición de redes hidráulicas y la resolución de las mismas. El estudio y el manejo de esta librería ha sido uno de los puntos claves del presente trabajo. Con ella se consigue generar la red de Alperovits y resolverla tantas veces como es necesario. Con el objetivo de comprobar si el algoritmo funciona correctamente y permite obtener buenas soluciones en relación a las alcanzadas anteriormente por otros autores, se simulan diferentes escenarios. La mejor solución encontrada tiene un coste de 420,000 u.m. A pesar de lograr buenos resultados en tiempos de resolución bajos, no se consigue alcanzar la solución de menor coste encontrada en la literatura científica, 419,000 u.m., 75 utilizando otros métodos de resolución. Además, no hay ninguna prueba de que esta solución sea un mínimo global del problema, es decir, es posible que dicha solución aun pueda optimizarse. Sería de especial interés la aplicación de la técnica diseñada a redes más complejas como las presentadas en este trabajo, red de Hanoi y de New York. Para comprobar si los tiempos de resolución aumentan y si las soluciones encontrar también son de calidad en relación a las encontradas por otros autores. Por último, se plantea como línea de trabajo a seguir en el futuro:  La modificación el algoritmo de búsqueda tabú añadiendo memoria a largo plazo.  El contraste de la solución alcanzada añadiendo términos de penalización a la función de evaluación de costes.  La generación con el paquete WNTR de las redes de Hanoi y de New York con el fin de verificar el funcionamiento del algoritmo de búsqueda Tabú diseñado y comparar los resultados con los obtenidos mediante otros métodos de resolución. 76 8. Referencias *1+ Datos del agua, Fundación Aquae. Enlace: https://www.fundacionaquae.org. Fecha de última consulta: 26.06.2018. *2+ Programa de las Naciones Unidas para el Medio Ambiente, titulado GEO4: Perspectivas del Medio Ambiente Mundial, 2006. Datos: Shikanov y Radoo, 2003. *3+ Instituto Español de Estadística, Agricultura y Medio Ambiente, Agua, Estadísticas sobre el suministro y saneamiento del agua, 2014, Código IOE 30077. *4+ Enlace: https://www.romanorumvita.com/. Fecha de última consulta: 06.07.2018. *5+ Instrucciones técnicas para redes de abastecimiento, PD 005 02, Revisión Nº4, EMASESA, 21 de Octubre de 2013. *6+ Jaime Navarro & Ana Prieto, Apuntes Redes de abastecimiento de agua en poblaciones, Asignatura: Servicios Urbanos, Plan 2011, E.T.S. Arquitectura, Universidad de Sevilla. *7+ E. Alperovits & U. Shamir, Design of Optimal Water Distribution Systems, Water Resources Research, 1977, Vol 13. No 6. *8+ O. Fujiwara & D.B. Khang, A Two-Phase Decomposition Method for Optimal Design of Looped Water Distribution Networks, Water Resources Research, 1990, Vol 26, Issue 4, pp. 539-549. *9+ Quindry, G . E., E. D. Brill, and J. C. Liebman, Optimization of looped water distribution systems, J. Environ. Eng. Div. Am. Soc. Civ. Eng., 1981, 107(EE4), pp. 665-679. *10+ C. R. Suribabu & T. R. Neelakantan, Design of water distribution networks using particle swarm optimization, Urban Water Journal, June 2006, Vol. 3, No. 2, pp. 111 – 120. *11+ Cristóbal Madero y de Miguel, Diseño y planificación de redes hidráulicas de distribución de Agua Potable, Trabajo fin de Máster, 2011. *12+ L.F. Moody, Friction Factors for Pipe Flow, Transactions of the ASME, 1944, vol 66, pp. 671-684. 77 *13+ Kessler, A., Shamir, U., Analysis of the linear programming gradient method for optimal design of water supply networks. Water Resources Research, 1989, 25 (7), pp. 1469–1480. *14+ Varma, K.V.K., Narasimhan, S., Bhallamudi, S.M., Optimal design of water distribution systems using an NLP method. Journal of Environmental Engineering, 1997, ASCE 123 (4), pp. 381–388. *15+ Cunha, M.C., Sousa, J., Water distribution network design optimization: Simulated annealing approach. Journal of Water Resources Planning and Management, 1999, ASCE 125 (4), pp. 215–221. *16+ Cunha, M.C., Sousa, J., Hydraulic infrastructures design using simulated annealing. Journal of Infrastructure Systems, 2001, ASCE 7 (1), pp. 32–39. *17+ Meirelles G., Melo B., Luvizotto Jr. E., Optimal design of water supply networks using an energy recovery approach, Renewable Energy 117, 2018, pp. 404-413. *18+ Matías, A., Diseño de redes de distribución de agua contemplando la fiabilidad mediante algoritmos genéticos, Ph.D. thesis, Politech, 2003, Univ. of Valencia, Valencia, Spain. *19+ Afshar, M.H., Application of a mix–min ant system to joint layout and size optimization of pipe networks. Engineering optimization, 2006, 38 (3), pp. 299 – 317. *20+ M. Cunha & L. Ribeiro, Tabu search algorithms for water network optimization, European Journal of Operational Research 157, 2004, pp. 746–758. *21+ Python Software Foundation. Downloaded from: https://www.python.org/ *22+ Fred Glover, Future paths for integer programming and links to artificial intelligence, Comput. & Ops. Res. Vol 13, No. 5, 1986, pp. 533-549. *23+ Alejandro Escudero & Pablo Aparicio, Apuntes de Modelado y optimización de problemas de gestión, E.T.S.I., Universidad de Sevilla, 2018. *24+ PyCharm, Python IDE for professional developers. Downloaded from: https://www.jetbrains.com/pycharm/ *25+ Klise, K.A., Bynum, M., Moriarty, D., Murray, R. , A software framework for assessing the resilience of drinking water systems to disasters with an example earthquake case study, Environmental Modelling y Software, 2017, Vol 95, pp. 420-431. 78 *26+ Klise, K.A., Hart, D.B., Moriarty, D., Bynum, M., Murray, R., Burkhardt, J., Haxton, T., Water Network Tool for Resilience (WNTR) User Manual, U.S. Environmental Protection Agency Technical Report, 2017, EPA/600/R-17/264, 47p. *27+ Rossman, L.A., EPANET, User’s Manual, 2000 (US Environmental Protection Agency: Cincinnati). *28+ R. Baños, C. Gil, J. Reca, F.G. Montoya, A memetic algorithm applied to the design of water distribution networks, Applied Soft Computing, 2010, Vol 10, pp. 261–266. *29+ Savic, D.A., Walters, G.A., Genetic algorithms for least cost design of water distribution networks. Journal of Water Resources Planning and Management, 1997, ASCE 123 (2), pp. 67–77. 79 ANEXO # Autor: Alicia Robles Velasco # Date: 03.07.2018 import wntr import pandas as pd import matplotlib.pylab as plt import numpy as np from time import time import random import math import operator import copy import pickle def Alperovits_Network(D): # --- Units --- # # demand: m3/s # coordinates: m # elevation: m # length: m # diameters: m Demand=[] Demand_in=[100,100,120,270,330,200] # In m3/h for i in range(len(Demand_in)): Demand.append(Demand_in[i]/3600) # Converted into m3/s # --- Create Network --- # wn = wntr.network.WaterNetworkModel() # --- Nodes --- # wn.add_reservoir('res', base_head=210, head_pattern=None, coordinates=(2000,2000)) wn.add_junction('node2', base_demand=Demand[0],elevation=150.0, coordinates=(1000,2000)) wn.add_junction('node3', base_demand=Demand[1],elevation=160.0, coordinates=(0,2000)) wn.add_junction('node4', base_demand=Demand[2],elevation=155.0, coordinates=(1000,1000)) wn.add_junction('node5', base_demand=Demand[3],elevation=150.0, coordinates=(0,1000)) wn.add_junction('node6', base_demand=Demand[4],elevation=165.0, coordinates=(1000,0)) wn.add_junction('node7', base_demand=Demand[5],elevation=160.0, coordinates=(0,0)) # --- Links --- # wn.add_pipe('pipe1', 'res', 'node2', length=1000, diameter=D[0], roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe2', 'node2', 'node3', length=1000, diameter=D[1], roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe3', 'node2', 'node4', length=1000, diameter=D[2], 80 roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe4', 'node4', 'node5', length=1000, diameter=D[3], roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe5', 'node4', 'node6', length=1000, diameter=D[4], roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe6', 'node6', 'node7', length=1000, diameter=D[5], roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe7', 'node3', 'node5', length=1000, diameter=D[6], roughness=130,minor_loss=0.0, status='OPEN') wn.add_pipe('pipe8', 'node5', 'node7', length=1000, diameter=D[7], roughness=130,minor_loss=0.0, status='OPEN') return wn def Decrease_diameter(D,i): # --- Commercial diameters [m] --- # Commercial_D=[0.0254,0.0508,0.0762,0.1016,0.1524,0.2032,0.254,0.30 48,0.3556,0.4064,0.4572,0.508,0.5588,0.6096] new_D=[] for j in range(len(D)): new_D.append(D[j]) if j == i: del new_D[j] if D[j]==Commercial_D[0]: new_D.append(D[j]) else: aux=[Commercial_D[z-1] for z in range(1,len(Commercial_D)) if D[j] == Commercial_D[z]] new_D.append(aux[0]) return new_D def Increase_diameter(D,i): # --- Commercial diameters [m] --- # Commercial_D=[0.0254,0.0508,0.0762,0.1016,0.1524,0.2032,0.254,0.30 48,0.3556,0.4064,0.4572,0.508,0.5588,0.6096] new_D=[] for j in range(len(D)): new_D.append(D[j]) if j == i: del new_D[j] if D[j]==Commercial_D[13]: new_D.append(D[j]) else: aux=[Commercial_D[z+1] for z in range(len(Commercial_D)-1) if D[j] == Commercial_D[z]] new_D.append(aux[0]) return new_D def Change_Network_Diameters(wn,D): for i in range(1,len(wn.links)+1): pipe=wn.get_link('pipe'+str(i)) pipe.diameter=D[i-1] return wn def Eval_Network_Technical(wn,D): # ----- Evaluation through WNTR package: EpanetSimulator ----- # min_H=30 wn=Change_Network_Diameters(wn,D) epanet_sim = wntr.sim.EpanetSimulator(wn) epanet_sim_results = epanet_sim.run_sim() Cont=0 81 # ----- Constrains: Pressure in the nodes ----- # for i in range (len(epanet_sim_results.node['pressure'].values[0])): if epanet_sim_results.node['pressure'].values[0][i]>min_H: Cont+=1 if Cont>=(len(epanet_sim_results.node['pressure'].values[0]) -1): Technical_Criterium=True else: Technical_Criterium=False return Technical_Criterium def Invalid(Invalid_solution,D): Invalid=False for i in range(len(Invalid_solution)): if D==Invalid_solution[i]: Invalid=True return Invalid def Eval_Network_Cost(D): # ----- Commercial diameters [m] ----- # Commercial_D=[0.0254,0.0508,0.0762,0.1016,0.1524,0.2032, 0.254,0.3048,0.3556,0.4064,0.4572,0.508,0.5588,0.6096] # ----- Cost of diameters [u.m./m] ----- # Cost_D=[2,5,8,11,16,23,32,50,60,90,130,170,300,550] Cost=0 # ----- Lenth of the arcs ----- # L=[1000,1000,1000,1000,1000,1000,1000,1000] for i in range(len(D)): for j in range(len(Commercial_D)): if D[i]==Commercial_D[j]: Cost+=(Cost_D[j]*L[i]) return Cost def Tabu_Search(D): # ----- Parameters ----- # tabu_size = len(D)//4 n_iteration=1000 Graph_BS=[]#To stock Best Solution Graph_CS=[]#To stock Current Solution Invalid_solutions=[] #To stock Invalid Solutions wn=Alperovits_Network(D) # ----- Initialitation ----- # current_solution=D[:] current_solution_value=Eval_Network_Cost(current_solution) best_solution=current_solution[:] best_solution_value=current_solution_value Graph_BS.append(best_solution_value) Graph_CS.append(current_solution_value) tabu_list=np.zeros([len(current_solution)],dtype=int) Invalid_solutions=[[0 for i in range(len(current_solution))]] for _ in range(n_iteration): # ----- Neighbourhood ----- # neighbourhood=[] neighbourhood_value=[] sec=[] for i in range(len(current_solution)): neighbour=Decrease_diameter(current_solution,i) 82 if Invalid(Invalid_solutions,neighbour)==False: if Eval_Network_Technical(wn,neighbour)==True: # ----- Aspiration Criterion ----- # if tabu_list[i]==0 or Eval_Network_Cost (neighbour)<best_solution_value: neighbourhood.append(neighbour) neighbourhood_value.append(Eval_Network _Cost(neighbour)) sec.append(i) else: Invalid_solutions.append(neighbour) if len(neighbourhood)>0: current_solution_value,sec,current_solution=min(zip (neighbourhood_value,sec,neighbourhood)) Graph_CS.append(current_solution_value) # ----- Update best solution ----- # if current_solution_value < best_solution_value: best_solution=current_solution[:] best_solution_value=current_solution_value Graph_BS.append(best_solution_value) # ----- Update tabu list ----- # tabu_list[sec]=tabu_size+1 for i in range(len(tabu_list)): if tabu_list[i]>0: tabu_list[i]-=1 else: # --- Any valid solution : Increase a diameter randomly --- # aux=[i for i in range(len(current_solution)) if tabu_list[i]==0] random.shuffle(aux) #Random choice current_solution=Increase_diameter( current_solution,aux[0]) current_solution_value=Eval_Network_Cost(current_solution) Graph_CS.append(current_solution_value) sec=aux[0] tabu_list[sec]=tabu_size+1 for i in range(len(tabu_list)): if tabu_list[i]>0: tabu_list[i]-=1 # ----- Graph of solutions ----- # Graph_solution(Graph_BS,Graph_CS) wn=Alperovits_Network(best_solution) Graph_network(wn) # ----- Velocity ----- # epanet_sim = wntr.sim.EpanetSimulator(wn) epanet_sim_results = epanet_sim.run_sim() print('Pressure results – Best_ Solution:',epanet_sim_results.node['pressure']) print('Flowrate results – Best Solution:',epanet_sim_results.link['flowrate']) velocity=[] velocity=[4*epanet_sim_results.link['flowrate'].values[0][i]/ ((best_solution[i]**2)*3.14) for i in range(len(best_solution))] 83 print('Velocity results - Best Solution:',velocity) # ---- Invalid Solutions ----- # print("Número de soluciones no validas encontradas:" ,len(Invalid_solutions)) return best_solution_value, best_solution def Graph_solution(Graph_BS,Graph_CS): colors = ['b', 'c', 'y', 'm', 'r', 'g', 'k', 'w'] plt.plot(Graph_CS,color=colors[4],linewidth=1) plt.title("Current_Solution_Value") plt.savefig("C:/Users/alici/Documents/Universidad/Máster/TFM/ Current_Solution.png") plt.show() plt.plot(Graph_BS,color=colors[4],linewidth=1) plt.title("Best_Solution_Value") plt.savefig("C:/Users/alici/Documents/Universidad/Máster/TFM/ Best_Solution.png") plt.show() def Graph_network(wn): epanet_sim = wntr.sim.EpanetSimulator(wn) epanet_sim_results = epanet_sim.run_sim() wntr.graphics.plot_network(wn, node_attribute = 'elevation', node_size=50, title='Node elevation') plt.savefig("C:/Users/alicia/Documents/Universidad/Máster/TFM/ Node elevation.png") wntr.graphics.plot_network(wn, node_attribute= epanet_sim_results.node['demand'].loc[0,:], node_size=50, title='Node demand') plt.savefig("C:/Users/alicia/Documents/Universidad/Máster/TFM/ Node demand.png") wntr.graphics.plot_network(wn, node_attribute= epanet_sim_results.node['pressure'].loc[0,:], node_size=50, title='Node pressure') plt.savefig("C:/Users/alicia/Documents/Universidad/Máster/TFM/ Node pressure.png") wntr.graphics.plot_network(wn, node_attribute = epanet_sim_results.node['head'].loc[0,:], node_size=50, title='Node head') plt.savefig("C:/Users/alicia/Documents/Universidad/Máster/TFM/ Node head.png") wntr.graphics.plot_network(wn, link_attribute = epanet_sim_results.link['flowrate'].loc[0,:], link_width=4, title='Link flowrate') plt.savefig("C:/Users/alicia/Documents/Universidad/Máster/TFM/ Link flowrate.png") plt.show() def main(): # ----- Initial solution A ----- # #D=[0.6096,0.6096,0.6096,0.6096,0.6096,0.6096,0.6096,0.6096] # ----- Initial solution B ----- # D=[0.508,0.4064,0.4064,0.254,0.0508,0.3048,0.3048,0.1016]