Planificación de movimiento de vehículos acuáticos de superficie basado en Optimización Bayesiana

Equation Chapter 1 Section 1 Proyecto Fin de Grado Ingeniería Electrónica, Robótica y Mecatrónica Planificación de movimiento de vehículos acuáticos de superficie basado en Optimización Bayesiana Autor: Pedro José Díaz García Tutor: Daniel Gutiérrez Reina Dpto. Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 iii Proyecto Fin de Grado Ingeniería Electrónica, Robótica y Mecatrónica Planificación de movimiento de vehículos acuáticos de superficie basado en Optimización Bayesiana Autor: Pedro José Díaz García Tutor: Daniel Gutiérrez Reina Profesor Contratado Dpto. Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 v Proyecto Fin de Grado: Planificación de movimiento de vehículos acuáticos de superficie basado en Optimización Bayesiana Autor: Pedro José Díaz García Tutor: Daniel Gutiérrez Reina El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2020 El Secretario del Tribunal vii A todos los que me han acompañado en este camino. ix Agradecimientos Debo empezar por el principio, y este trabajo comienza con todas aquellas personas que me han brindado su apoyo, a las que quiero expresar mi más sincero agradecimiento: A mi tutor Daniel Gutiérrez y a Federico Peralta por toda la ayuda prestada durante el desarrollo de este trabajo. A mi familia por ser un punto de apoyo indiscutible tanto en el ámbito académico como en el personal. Gracias por creer en mí siempre. A todos mis amigos, tanto fuera como dentro de clase, por soportarme en mis malos ratos, animarme en los peores y acompañarme en los mejores. En especial a Raúl, junto a quien volvería a recorrer este camino mil veces más. Una vez más, porque nunca es suficiente, gracias a todos. Pedro José Díaz García Sevilla, 2020 ÍNDICE DE TABLAS Tabla 4-1. Resumen de resultados para el experimento 1 28 Tabla 4-2. Resumen de resultados para el experimento 2 30 Tabla 4-3. Resumen de resultados para el experimento 3 32 Tabla 4-4. Resumen de resultados para el experimento 4 34 Tabla 4-5. Resumen de resultados para el experimento 5 36 Tabla 4-6. Resumen de resultados para el experimento 6 38 Tabla 4-7. Resumen de resultados para el experimento 7 40 Tabla 4-8. Resumen de resultados para el experimento 8 42 Tabla 4-9. Resumen de resultados para el experimento 9 44 Tabla 4-10. Resumen de resultados para el experimento 10 46 Tabla 4-11. Resumen de resultados para el experimento 11 48 Tabla 4-12. Resumen de resultados para el experimento 12 50 Tabla 4-13. Resumen de resultados para el experimento 13 52 Tabla 4-14. Resumen de resultados para el experimento 14 54 Tabla 4-15. Resumen de resultados para el experimento 15 56 Tabla 4-16. Resumen de resultados para el experimento 16 58 Tabla 4-17. Resumen de resultados para el experimento 17 60 Tabla 4-18. Resumen de resultados para el experimento 18 62 xvii ÍNDICE DE FIGURAS Figura 1-1. Geografía del lago Ypacaraí. A la izquierda, mapa de Paraguay, donde se indican las localizaciones del lago Ypacaraí y de la capital Asunción. A la derecha, mapa físico del lago Ypacaraí. Fuentes: https://upload.wikimedia.org/wikipedia/commons/e/e9/Paraguay_rel_location_map.svg 1 Figura 1-2. Floración de cianobacterias. 2 Figura 1-3. Prototipos de ASV. A la izquierda, Cormorán II. A la derecha, Cormorán III. Fuente: [3] 3 Figura 2-1. Representación de la posición de las balizas sobre el contorno del lago Ypacaraí. Fuente: [10] 7 Figura 3-1. Representación gráfica para el caso de una dimensión del funcionamiento de la Optimización Bayesiana. Fuente: [14] 12 Figura 3-2. Modelado mediante procesos gaussianos. Fuente: [32] 14 Figura 3-3. Representación de proceso gaussiano. En cada punto, el proceso gaussiano devuelve los valores de media y desviación típica que definen la distribución normal que siguen los valores de la función en dicho punto. Fuente: [14] 14 Figura 3-4. Variación del kernel RBF en función del valor del lenght-scale. Arriba l=1.5, centro l=1, abajo l=0.5. Fuente: https://distill.pub/2019/visual-exploration-gaussian-processes/ 16 Figura 4-1. Representación tridimensional de la función de Branin. Fuente: [27] 21 Figura 4-2. Representación tridimensional de la función de Schwefel (para d=2). Fuente: [27] 22 Figura 4-3. Representación tridimensional de la función de Styblinski-Tang (para d=2). Fuente: [27] 23 Figura 4-4. Diagrama de flujo del código implementado. 25 Figura 4-5. Resultados obtenidos para el experimento 1 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 28 Figura 4-6. Evolución de los mapas de error y precisión durante el experimento 1 29 Figura 4-7. Resultados obtenidos para el experimento 2 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 30 Figura 4-8. Evolución de los mapas de error y precisión durante el experimento 2 31 Figura 4-9. Resultados obtenidos para el experimento 3 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 32 Figura 4-10. Evolución de los mapas de error y precisión durante el experimento 3 33 Figura 4-11. Resultados obtenidos para el experimento 4 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 34 Figura 4-12. Evolución de los mapas de error y precisión durante el experimento 4 35 Figura 4-13. Resultados obtenidos para el experimento 5 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 36 Figura 4-14. Evolución de los mapas de error y precisión durante el experimento 5 37 Figura 4-15. Resultados obtenidos para el experimento 6 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 38 xix Figura 4-16. Evolución de los mapas de error y precisión durante el experimento 6 39 Figura 4-17. Resultados obtenidos para el experimento 7 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 40 Figura 4-18. Evolución de los mapas de error y precisión durante el experimento 7 41 Figura 4-19. Resultados obtenidos para el experimento 8 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 42 Figura 4-20. Evolución de los mapas de error y precisión durante el experimento 8 43 Figura 4-21. Resultados obtenidos para el experimento 9 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 44 Figura 4-22. Evolución de los mapas de error y precisión durante el experimento 9 45 Figura 4-23. Resultados obtenidos para el experimento 10 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 46 Figura 4-24. Evolución de los mapas de error y precisión durante el experimento 10 47 Figura 4-25. Resultados obtenidos para el experimento 11 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 48 Figura 4-26. Evolución de los mapas de error y precisión durante el experimento 11 49 Figura 4-27. Resultados obtenidos para el experimento 12 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 50 Figura 4-28. Evolución de los mapas de error y precisión durante el experimento 12 51 Figura 4-29. Resultados obtenidos para el experimento 13 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 52 Figura 4-30. Evolución de los mapas de error y precisión durante el experimento 13 53 Figura 4-31. Resultados obtenidos para el experimento 14 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 54 Figura 4-32. Evolución de los mapas de error y precisión durante el experimento 14 55 Figura 4-33. Resultados obtenidos para el experimento 15 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 56 Figura 4-34. Evolución de los mapas de error y precisión durante el experimento 15 57 Figura 4-35. Resultados obtenidos para el experimento 16 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 58 Figura 4-36. Evolución de los mapas de error y precisión durante el experimento 16 59 Figura 4-37. Resultados obtenidos para el experimento 17 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 60 Figura 4-38. Evolución de los mapas de error y precisión durante el experimento 17 61 Figura 4-39. Resultados obtenidos para el experimento 18 a las n iteraciones. De izquierda a derecha: ξ=0.01, ξ=0.1, ξ=0.2 62 Figura 4-40. Evolución de los mapas de error y precisión durante el experimento 18 63 Figura 4-41. Resumen de los resultados de error para n iteraciones. 64 Figura 5-1. Esquema del sistema GNC del ASV. Fuente: [28] 67 xxi Notación km Kilómetros km2 Kilómetros cuadrados m Metros p. ej. Por ejemplo ∝ Relación de proporcionalidad exp Función Exponencial 𝜋 Número Pi GP Proceso gaussiano (Gaussian Process) < Menor que > Mayor que 𝒩 Distribución Normal argmax Argumento del máximo max Máximo 1 INTRODUCCIÓN 1.1 Antecedentes El lago Ypacaraí constituye una de las zonas geográficas más importantes de Paraguay. Ya sea por su atractivo turístico, o por el uso de sus aguas para actividades agrícola-ganaderas e industriales, el lago es una pieza clave en el motor económico del país. Dada la importancia que tiene, la calidad de sus aguas no desmerece una atención especial. Esta masa de agua dulce se sitúa a unos 37 km al este de la capital Asunción. Con un área total de unos 90 km2 y una longitud de 24 km de norte a sur y 5 a 6 km de este a oeste, se posiciona como el lago más grande de Paraguay. Se trata de un lago poco profundo, teniendo una profundidad media de 3 m [1]. El lago Ypacaraí cumple un papel clave para la sociedad de la zona. Por una parte, funciona como fuente de abastecimiento de agua para los municipios cercanos. Las distintas poblaciones de la zona aprovechan esta agua tanto para consumo directo como para llevar a cabo actividades ganaderas, agrícolas e industriales. Por otro lado, también constituye un importante enclave turístico con especial relevancia en la temporada estival, durante la cual sus aguas acogen tanto a bañistas como a practicantes de diferentes deportes acuáticos. Lago Ypacaraí Asunción Figura 1-1. Geografía del lago Ypacaraí. A la izquierda, mapa de Paraguay, donde se indican las localizaciones del lago Ypacaraí y de la capital Asunción. A la derecha, mapa físico del lago Ypacaraí. Fuentes: https://upload.wikimedia.org/wikipedia/commons/e/e9/Paraguay_rel_location_map.svg Introducción 2 2 En las últimas décadas, la presencia del lago ha favorecido el desarrollo económico y demográfico de la región. Sin embargo, esto ha conllevado un incremento de los residuos vertidos al lago a través de sus arroyos: desechos domésticos, urbanos e industriales, así como otros aportes difusos, en su mayoría inorgánicos, procedentes de la actividad agrícola-ganadera. Como consecuencia, ha tenido lugar un fenómeno conocido como eutrofización [2], en el que el vertido incontrolado de estos desechos da lugar a una acumulación excesiva de nutrientes, principalmente fósforo y nitrógeno, impidiendo que se lleve a cabo un ciclo de asimilación natural de estos. Adicionalmente, se produce una caída drástica en los niveles de oxígeno, provocando una gran disminución de la calidad del agua, haciéndola no apta para algunos de sus usos actuales, como son su uso doméstico o como escenario para actividades acuáticas. Adicionalmente, la presencia de estos nutrientes en altas concentraciones favorece la aparición de cianobacterias, dando a las aguas un color verde frente al azulado que tenían anteriormente. El principal peligro de la aparición de cianobacterias es que la combinación de las floraciones de este tipo de alga junto a las cianotoxinas que producen puede ser un problema tanto para el ser humano como para otros organismos que habiten el medio acuático. Aunque los niveles detectados de cianobacterias están por debajo de los niveles mínimos para que sea considerado peligroso, estos niveles deben ser tratados ya que estas algas producen un fuerte y desagradable olor que perjudica la imagen del lago, resultando poco atractivo para el ocio de sus usuarios. Con la finalidad de solventar este problema, se han llevado a cabo diferentes proyectos para tomar medida de la calidad del agua. Sin embargo, estos resultaban ser poco eficientes, ya que las mediciones eran poco continuas y se realizaban en puntos específicos, además de que requerían de la presencia de personal técnico para tomar las muestras. Es en 2016 cuando arranca el presente proyecto, con la promesa de innovar en los aspectos en los que los anteriores proyectos fallaban. Este trabajo se enmarca en la línea de cooperación entre la Universidad de Sevilla (US) y la Universidad Nacional de Asunción (UNA) para la monitorización de la calidad del agua en el lago Ypacaraí. En el caso de este proyecto, se pretende llevar a cabo la toma de datos empleando una flota de drones acuáticos de superficie que recorrerá la mayor superficie del lago posible recabando información en diferentes puntos de esta. De esta manera, no será necesario personal técnico presencial para tomar las muestras ni instalaciones estáticas. Figura 1-2. Floración de cianobacterias. Fuente: https://www.hoy.com.py/nacionales/vuelve-la-espuma-verde-al-lago-ypacarai 3 En concreto, el proyecto de colaboración entre la UNA y la US persigue cuatro objetivos fundamentales: - Construir una red de vehículos no tripulados provistos de sensores que les permitan captar información útil para medir la calidad del agua, y con capacidad de comunicación entre ellos mediante una red inalámbrica que les permita compartir la información capturada. - Otorgar a la red de vehículos la capacidad de aprendizaje continuo sobre el entorno cambiante, actualizando los modelos a medida que se incorporan nuevos datos. - Poner en marcha un sistema de telemetría y monitorización en tiempo real de las condiciones del lago que permita determinar cómo influye el entorno sobre las floraciones de cianobacterias en base a los datos capturados. - Poner el sistema desarrollado a disposición de las entidades locales encargadas de la gestión del agua del lago. La monitorización continua del lago permitirá mejorar la información disponible para la población, la sensibilización sobre el problema y un posible desarrollo normativo que lo solucione o mitigue. El proyecto en todo su conjunto resulta innovador en cuanto a la medición de contaminación en ríos y lagos, pero también ofrece la posibilidad de ser aplicado a otros problemas ambientales. No solo eso, sino que las investigaciones aquí realizadas serán extrapolables a otras aplicaciones relacionadas con el muestreo mediante vehículos autónomos no tripulados. En el caso de este proyecto, se está trabajando con dos prototipos de vehículo autónomo de superficie (ASV, Autonomous Surface Vehicle): el Cormorán II y el Cormorán III. En el caso del primero, se trata de un vehículo acuático de 4x2 m, provisto de dos sistemas de detección de obstáculos (un LIDAR y un sistema de procesamiento de imágenes basado en redes neuronales convolucionales), y alimentado mediante dos baterías de 20 Ah y paneles fotovoltaicos que le permiten alargar su autonomía. En el caso del Cormorán III, se trata de un vehículo de menores dimensiones (1.3 x 0.95 m), cuyos sensores y electrónica de control son los mismos que se tenían en el Comorán II, aunque se pierden algunas prestaciones, como el uso de placas fotovoltaicas, lo que reduce su autonomía. En caso de que el lector quiera saber más acerca de los ASV desarrollados para este proyecto, puede consultar [3]. Figura 1-3. Prototipos de ASV. A la izquierda, Cormorán II. A la derecha, Cormorán III. Fuente: [3] Optimización Bayesiana 10 10 3 OPTIMIZACIÓN BAYESIANA La Optimización Bayesiana es un método de optimización basado en Machine Learning cuya principal aplicación tiene lugar sobre funciones objetivo con alguna de las siguientes características: - La función objetivo es a priori desconocida, usualmente denominada black-box function, lo que conlleva varios inconvenientes. Por una parte, imposibilita conocer las derivadas de la función, impidiendo el uso de métodos de optimización clásicos como el descenso del gradiente o el método de Newton. Por otro lado, limita las suposiciones que pueden hacerse sobre la función y que ayudarían a mejorar la eficiencia del proceso de optimización, como su linealidad o su convexidad. De hecho, al no conocer su convexidad, debe presuponerse la existencia de máximos o mínimos locales, forzando a utilizar un método de optimización global, con la consecuente delimitación del espacio a optimizar. - La función objetivo es costosa de evaluar en el sentido de que el número de observaciones que se pueden realizar sobre esta no es ilimitado. Esta limitación puede venir dada porque el coste monetario de realizar la observación sea elevado (p.ej.: ensayos destructivos), porque el tiempo que conlleve realizar la observación sea demasiado elevado (p.ej.: simulaciones complejas), o porque las oportunidades de realizar observaciones sean limitadas. - La función objetivo puede presentar ruido. Si bien es cierto que la Optimización Bayesiana es perfectamente aplicable a casos libres de ruido, esta técnica se presenta como una alternativa robusta en casos donde el ruido sí está presente. Según esto, el objetivo de la Optimización Bayesiana es encontrar aquellos valores de la función que aporten la mayor cantidad de información posible sobre esta. Esto permitirá que el número de observaciones necesarias sea menor, de forma que se reduce el tiempo y coste de optimización. Es importante tener en cuenta las condiciones bajo las que se debe hacer uso de esta técnica si se quieren obtener buenos resultados, ya que será en estos casos donde la Optimización Bayesiana alcance su máximo potencial. En comparación con otras técnicas de optimización, como pueden ser los algoritmos genéticos o la optimización por enjambre, la Optimización Bayesiana, por lo general, obtiene peores resultados. Sin embargo, en el caso de que la función objetivo cumpla alguna de las condiciones expuestas anteriormente, dichas técnicas alternativas resultan muy ineficientes o imposibles de aplicar, ya que basan su funcionamiento en evaluar muchas veces la función objetivo, siendo esto en ocasiones muy costoso de realizar. 11 3.1 Fundamento y algoritmo Tal como su nombre indica, la Optimización Bayesiana tiene su base en el Teorema de Bayes (Thomas Bayes, 1763), que relaciona las probabilidades de que ocurra un evento aleatorio A condicionado a otro evento aleatorio B y la de que ocurra B condicionado a A. El teorema se define como sigue [26]: Teorema 2-1 (Teorema de Bayes) Sea {𝐴1,𝐴2,…,𝐴𝑖,…,𝐴𝑛} un conjunto de sucesos mutuamente excluyentes y exhaustivos, y tales que la probabilidad de cada uno de ellos es distinta de cero. Sea B un suceso cualquiera del que se conocen las probabilidades condicionales 𝑃(𝐵|𝐴𝑖). Entonces, la probabilidad 𝑃(𝐴𝑖|𝐵) viene dada por la expresión: 𝑃(𝐴𝑖|𝐵)=𝑃(𝐵|𝐴𝑖)𝑃(𝐴𝑖) 𝑃(𝐵) donde: 𝑃(𝐴𝑖) y 𝑃(𝐵) son las probabilidades a priori de los sucesos 𝐴𝑖 y 𝐵, respectivamente. 𝑃(𝐵|𝐴𝑖) es la probabilidad de B condicionada a 𝐴𝑖. 𝑃(𝐴𝑖|𝐵) es la probabilidad a posteriori. Aplicado el teorema anterior al caso de un problema de optimización, la expresión anterior toma la forma: 𝑃(𝑀|𝐸)∝𝑃(𝐸|𝑀)𝑃(𝑀) (3-1) donde M denota al modelo y E a la evidencia, es decir, al conjunto de observaciones realizadas. Nótese que, en comparación con la definición formal del teorema, ha desaparecido el término del denominador, por lo que deben entenderse las nuevas probabilidades como “probabilidades normalizadas” por dicho valor. Por ello, se pierde la condición de igualdad, pasando a reconocerse como una relación de proporcionalidad; aunque será igualmente válida para el objetivo perseguido. La expresión (3-1) puede interpretarse como: La probabilidad a posteriori del modelo M dada una evidencia E es proporcional al producto de la probabilidad de obtener dicha evidencia E dado M y la probabilidad a priori del modelo M. Para llevar a cabo un muestreo eficiente de la función de evaluación aplicando el teorema de Bayes, se sustituye la función objetivo por un modelo más sencillo (modelo subrogado), denominado prior. Este modelo otorgará información sobre la probabilidad de los datos, y se irá mejorando en base a las evidencias que se obtengan del modelo real en el proceso de muestreo, entendiéndose como mejora un mayor acercamiento a la función real que se está modelando. De esta forma, se reescribe la expresión anterior como: 𝑃(𝑓|𝐷1:𝑡)∝𝑃(𝐷1:𝑡|𝑓)𝑃(𝑓) (3-2) donde f es la función objetivo real y 𝐷1:𝑡={𝒙1:𝑡,𝑓(𝒙1:𝑡)} es el conjunto de muestras tomadas sobre dicha función. Como se ha especificado anteriormente, esta técnica se aplica cuando la función objetivo es costosa de evaluar. Por ello, para llevar a cabo la optimización, se utiliza otra función, denominada función de adquisición. Esta función estará basada en el modelo subrogado anterior, pero será más simple, conllevando un coste de evaluación mucho menor en comparación con la función objetivo real. Optimización Bayesiana 12 12 De esta forma, el proceso de optimización se realizará sobre la función de adquisición en lugar de sobre la función objetivo, y se tomará la siguiente muestra del modelo real en el punto donde la función de adquisición sea óptima. Con lo expuesto hasta ahora, puede desarrollarse el siguiente algoritmo para aplicar la técnica de Optimización Bayesiana: Algoritmo 2-1 (Optimización Bayesiana) for 𝑡=1,2,… do 1. Encontrar el siguiente punto a muestrear 𝒙𝑡 optimizando la función de adquisición 𝑢: 𝒙𝑡=𝑎𝑟𝑔𝑚𝑎𝑥𝑥 𝑢(𝒙|𝐷1:𝑡−1) 2. Muestrear la función objetivo en dicho punto: 𝑦𝑡=𝑓(𝒙𝑡)+𝜀𝑡 3. Actualizar el conjunto de datos 𝐷1:𝑡={𝐷1:𝑡−1;(𝒙𝑡,𝑦𝑡)} con las nuevas muestras. 4. Actualizar el modelo con los nuevos datos. end for Figura 3-1. Representación gráfica para el caso de una dimensión del funcionamiento de la Optimización Bayesiana. Fuente: [14] 13 Puede apreciarse que el algoritmo anterior se divide en dos fases: en primer lugar, se llevan a cabo la optimización y el muestreo y, posteriormente, se realiza la actualización del modelo con los nuevos datos. De esta forma, al final de cada iteración, se tiene el modelo actualizado. A estas fases se les denomina, respectivamente, fase de decisión y fase de aprendizaje. Empleando este algoritmo, en cada iteración el modelo subrogado irá mejorando, por lo que los valores óptimos obtenidos al optimizar la función de adquisición serán más parecidos a los valores óptimos de la función real. Obviamente, los resultados serán buenos en tanto que el modelo subrogado que se este utilizando sea adecuado para el problema. Al utilizar el modelo para guiar el proceso de optimización, la Optimización Bayesiana obtiene mejores resultados que otros métodos, como la búsqueda en rejilla o la búsqueda aleatoria. A modo de resumen, se han presentado dos herramientas clave en la implementación de la Optimización Bayesiana: el prior, que generalmente será modelado mediante un proceso gaussiano, y la función de adquisición; la obtención de buenos resultados pasará por una correcta elección de ambas herramientas. Adicionalmente, es típico encontrar hiperparámetros en la definición del prior por lo que será conveniente llevar a cabo un proceso previo de optimización para encontrar los valores óptimos de dichos hiperparámetros. 3.2 Procesos Gaussianos Como se introdujo anteriormente, la función objetivo real se sustituirá por un modelo más sencillo. Esta simplificación conlleva que la obtención de buenos resultados estará condicionada a la fidelidad con la que dicho modelo represente el comportamiento de la función real. Esta fidelidad sería mayor en tanto mayor fuera el número de propiedades asumibles de la función real, con el inconveniente de que esta puede ser desconocida. De cara a afrontar este problema de modelado, se dan dos puntos de vista comunes. Un primer acercamiento pasaría por restringir el tipo de funciones que se debe considerar (p.ej.: funciones lineales), descartando todos los demás casos. Esta metodología tiene un claro problema: en caso de escoger un tipo de función que no se asemeje a la función subyacente, las predicciones realizadas serán erróneas. El segundo acercamiento pretendería mitigar este inconveniente, siendo mucho más flexible: en lugar de restringir el tipo de función, se otorgaría una probabilidad a cada función, siendo esta mayor cuanto más semejante a la función subyacente se considere (p.ej.: otorgar mayor probabilidad a funciones con mayor suavidad). Sin embargo, esta solución también presenta un serio problema: cómo tener en consideración las infinitas funciones posibles y, más allá, cómo computarlas en un tiempo finito. A efectos de salvar estos inconvenientes, los procesos gaussianos constituyen una herramienta perfecta. Un proceso gaussiano es una extensión multivariante de la distribución gaussiana, definiéndose como una función de probabilidad sobre funciones aleatorias. Es decir, el proceso gaussiano puede interpretarse como una función que devuelve en cada punto una media y una varianza de una distribución normal, adquiriendo una gran flexibilidad a la hora de ajustarse a cualquier tipo de función. También solventa el problema de la tratabilidad computacional, ya que no es necesario operar con los infinitos puntos de la función para observar las propiedades de un determinado grupo de ellos. Además, al no tratarse de un modelo paramétrico, no hay que preocuparse por si el modelo puede ajustarse o no a los nuevos datos; aunque se añadan más datos al modelo, este será suficientemente flexible para adaptarse a ellos. Optimización Bayesiana 14 14 En la Figura 2-2a puede verse cómo ante la falta de datos sobre la función real son muchas las funciones que pueden corresponderse con esta. Sin embargo, en la Figura 2-2b, tras el muestreo de algunos puntos de la función real, el número de funciones que pueden modelarla será menor, ya que solo serán válidas aquellas funciones que se ajusten a los puntos observados. Matemáticamente, al igual que la distribución gaussiana o normal queda definida en base a su media y su desviación típica tal y como indica la expresión (3-3), un proceso gaussiano se define en base a una función de media 𝑚(𝒙) y una función de covarianza 𝑘(𝒙,𝒙′), como indica la expresión (3-4). 𝑁(𝑥:𝜇,𝜎)=1 √2𝜋𝜎2𝑒𝑥𝑝[−12(𝑥−𝜇)2 𝜎2] (3-3) 𝑓(𝑥)~𝐺𝑃(𝑚(𝒙),𝑘(𝒙,𝒙′)) (3-4) Figura 3-2. Modelado mediante procesos gaussianos. Fuente: [32] Figura 3-3. Representación de proceso gaussiano. En cada punto, el proceso gaussiano devuelve los valores de media y desviación típica que definen la distribución normal que siguen los valores de la función en dicho punto. Fuente: [14] 15 Por lo general, la función de media se considera nula, por lo que el proceso gaussiano queda enteramente definido a partir de la elección de la función de covarianza, comúnmente denominada kernel. Será el kernel el que determine las propiedades del proceso gaussiano, pudiendo controlar a partir de dicha función cómo variará el prior a partir de los puntos muestreados. De esta forma, los valores de un proceso gaussiano con estas características seguirán una distribución normal multivariante: 𝒩(0,𝑲),𝑠𝑖𝑒𝑛𝑑𝑜 𝑲=[𝑘(𝒙1,𝒙1) ⋯ 𝑘(𝒙1,𝒙𝑡) ⋮ ⋱ ⋮ 𝑘(𝒙𝑡,𝒙1) … 𝑘(𝒙𝑡,𝒙𝑡)] (3-5) donde debe tenerse en cuanta que 𝑘(𝒙𝑖,𝒙𝑖)=1. Típicamente, los kernels tienen la propiedad de que puntos cercanos en el espacio de entrada estarán más correlacionados, es decir: 𝑆𝑖 ‖𝒙1−𝒙2‖< ‖𝒙1−𝒙3‖, 𝑒𝑛𝑡𝑜𝑛𝑐𝑒𝑠 𝑘(𝒙1,𝒙2)> 𝑘(𝒙1,𝒙2) Dado que el prior debe permanecer actualizado, debe obtenerse una expresión que permita definir el proceso gaussiano en base al proceso gaussiano de un estado anterior y los nuevos datos obtenidos en cada iteración. Partiendo de la expresión (3-5), el proceso gaussiano puede actulizarse según: [𝒇𝑡 𝑓𝑡+1]=𝒩(𝟎,[𝑲 𝒌 𝒌𝑻𝑘(𝑥𝑡+1,𝑥𝑡+1)]) (3-6) donde la matriz K ya se definió anteriormente y el vector k se construye como: 𝒌=[𝑘(𝑥𝑡+1,𝑥1) … 𝑘(𝑥𝑡+1,𝑥𝑡)] Bajo estas definiciones, se tienen las siguientes expresiones mediante las cuales se puede obtener la distribución del siguiente punto a muestrear: 𝑃(𝑓𝑡+1|𝐷1:𝑡,𝒙𝑡+1)=𝒩(𝜇𝑡(𝒙𝑡+1),𝜎𝑡2(𝒙𝑡+1)) {𝜇𝑡(𝒙𝑡+1)=𝒌𝑇𝑲−1𝒇1:𝑡 𝜎𝑡2(𝒙𝑡+1)=𝑘(𝒙𝑡+1,𝒙𝑡+1)− 𝒌𝑇𝑲−1𝒌 (3-7) En la expresión anterior, resulta destacable que el cálculo queda completamente definido en base a la función de covarianza empleada, por lo que la elección de esta resulta clave para definir el proceso gaussiano. A continuación, se presentan dos de los kernels más utilizados. Optimización Bayesiana 16 16 Radial-basis function kernel (RBF) 𝑘(𝒙𝑖,𝒙𝑗)=𝑒𝑥𝑝(−1 2𝑙2‖𝒙𝑖−𝒙𝑗‖2) (𝑙>0) (3-8) El kernel RBF, también conocido como squared exponential kernel, es probablemente el kernel más usado en este campo, pese a que su alta suavidad resulta poco realista para modelar ciertos fenómenos. Esta función de covarianza queda parametrizada por 𝑙, denominado length-scale, que controla la flexibilidad del kernel a la hora de adaptarse a los datos: a menor 𝑙, mayor flexibilidad. Màtern kernel 𝑘(𝒙𝑖,𝒙𝑗)=1 2𝜈−1Γ(𝜈)(√2𝜈 𝑙‖𝒙𝑖−𝒙𝑗‖)𝜈𝐾𝜈(√2𝜈 𝑙‖𝒙𝑖−𝒙𝑗‖) (3-9) El grupo de los Màtern kernels consiste en una generalización del kernel RBF descrito anteriormente. En la expresión (3-9), Γ(∙) es la función gamma y 𝐾𝜈(∙) es una modificación de la función de Bessel. Por otra parte, este kernel cuenta con un parámetro adicional 𝜈 que permite controlar la suavidad de la función de covarianza resultante, siendo esta menos suave para valores de 𝜈 más pequeños. Figura 3-4. Variación del kernel RBF en función del valor del lenght-scale. Arriba l=1.5, centro l=1, abajo l=0.5. Fuente: https://distill.pub/2019/visual-exploration-gaussian-processes/ 17 De hecho, cuando 𝜈→∞ el kernel resultante es equivalente al kernel RBF, y cuando 𝜈=0.5 se convierte en el kernel exponencial absoluto. Además, existen otros valores típicos, como 𝜈=1.5 y 𝜈=2.5, para los que el costo computacional se ve reducido debido a que no es necesario reevaluar la función de Bessel para dichos valores. Los kernels anteriores, así como tantos otros que no han sido expuestos en este trabajo, cuentan con hiperparámetros (como 𝑙) cuyo valor deberá escogerse adecuadamente para un buen desempeño del algoritmo. Si bien siempre existe la posibilidad de escoger estos valores manualmente, también pueden aplicarse otros métodos que permitan determinar los valores óptimos para estos hiperparámetros y obtener así resultados mejores. Una característica excepcional de los procesos gaussianos es que, a diferencia de otros modelos de Machine Learning, cuentan con una expresión analítica computacionalmente tratable para la integral del modelo sobre los parámetros, lo que permite el cálculo del valor óptimo del hiperparámetro que permitirá el mejor ajuste del modelo. Esta expresión de denomina verosimilitud marginal, y viene dada por: 𝑝(𝒚|𝑋)=∫𝑝(𝒚|𝒇,𝑋)𝑝(𝒇|𝑋)𝑑𝒇 (3-10) En su forma logarítmica: log𝑝(𝒚|𝑋)=−12𝒚𝑻(𝐾+𝜎𝑛2𝐼)−1𝒚−12𝑙𝑜𝑔|𝐾+𝜎𝑛2𝐼|−𝑛2log (2𝜋) (3-11) Teniendo en cuenta que el kernel varía en base al hiperparámetro 𝜃, puede escribirse la siguiente expresión que permite llevar a cabo la optimización para obtener el valor óptimo del hiperparámetro: ∂ ∂θlog𝑝(𝒚|𝑋,𝜃)=12𝒚𝑻𝐾−1𝜕𝐾 𝜕𝜃𝐾−1𝒚−12𝑡𝑟(𝐾−1𝜕𝐾 𝜕𝜃) (3-12) 3.3 Funciones de adquisición Una vez se tiene el prior que modela la función objetivo real, se debe seleccionar una función de adquisición adecuada. Como se dijo anteriormente, será esta la función sobre la que se aplique el método de optimización, para hallar así el siguiente punto a muestrear. En efecto, la función de adquisición debe estar estrechamente relacionada con el prior, de manera que el punto óptimo encontrado se corresponda con un punto de alto interés en la función objetivo, ya sea por encontrarse cerca del óptimo real o por encontrarse en zonas inexploradas de la función objetivo. Es en este punto donde entra en juego el dilema exploración vs explotación. Se denomina exploración a la obtención de nuevas muestras en zonas de alta incertidumbre, es decir, en puntos no cercanos a muestras anteriores. Por su parte, se denomina explotación a la obtención de muestras en zonas donde se conoce que la función objetivo alcanza valores altos, es decir, en puntos cercanos a los puntos óptimos. Será la función de adquisición la que determine el tipo de tendencia que seguirá el algoritmo, por lo que la correcta elección de esta resulta de gran importancia para la obtención de los resultados perseguidos. Optimización Bayesiana 18 18 Con todo esto, la finalidad de la función de adquisición es doble: por una parte, debe permitir discernir cuáles son los mejores puntos a muestrear en el siguiente paso que permitan la mayor mejora del modelo y, por otro lado, debe estar definida por una expresión sencilla en cuanto a complejidad de cálculo, lo que permitirá el uso de métodos de optimización clásicos sobre las mismas. A continuación, se presentan algunas de las funciones de adquisición más utilizadas. Probability of Improvement (PI) Esta función mide la probabilidad de que la siguiente muestra sea mejor que el mejor valor obtenido hasta el momento. 𝑃𝐼(𝒙)=𝑃(𝑓(𝒙)≥𝑓(𝒙+))=Φ(𝜇(𝒙)−𝑓(𝒙+) 𝜎(𝒙) ),𝑐𝑜𝑛 𝒙+=𝑎𝑟𝑔𝑚𝑎𝑥𝑥𝑖∈𝑥1:𝑡 𝑓(𝒙𝑖) (3-13) donde 𝒙+ denota el punto donde se ha obtenido la mejor muestra hasta el momento, y Φ(∙) es la función de distribución acumulada normal. El mayor inconveniente de esta función es que tiende puramente a la explotación, es decir, una vez se encuentra un punto con un buen resultado, las siguientes muestras se tomarán en las cercanías de ese punto. En el caso de que se quiera que la mejora buscada supere un cierto umbral, se puede añadir un hiperparámetro que haga las veces de valor umbral: 𝑃𝐼(𝒙)=𝑃(𝑓(𝒙)≥𝑓(𝒙+)+𝜉)=Φ(𝜇(𝒙)−𝑓(𝒙+)−𝜉 𝜎(𝒙) ) (3-14) Expected Improvement (EI) En este caso, en lugar de utilizar la probabilidad de mejora, se utiliza la magnitud de la mejora esperada. 𝐸𝐼(𝒙)=Ε 𝑚𝑎𝑥(𝑓(𝒙)−𝑓(𝒙+),0) (3-15) Analíticamente, la mejora esperada puede ser evaluada como: 𝐸𝐼(𝒙)={(𝜇(𝒙)−𝑓(𝒙+)−𝜉)Φ(Z)+𝜎(𝒙)𝜙(𝑍) 𝑠𝑖 𝜎(𝒙)>0 0 𝑠𝑖 𝜎(𝒙)=0 (3-16) Donde Φ(∙) y 𝜙(∙) son la función de distribución acumulada normal (CDF, Cumulative Distribution Function) y la función de distribución de probabilidad normal (PDF, Probability Distribution Function), respectivamente. Por su parte, la distribución Z queda definida como: 𝑍={𝜇(𝒙)−𝑓(𝒙+)−𝜉 𝜎(𝒙) 𝑠𝑖 𝜎(𝒙)>0 0 𝑠𝑖 𝜎(𝒙)=0 19 En la expresión de la mejora esperada anterior pueden diferenciarse dos términos sumados. El primero se denomina término de explotación, y el segundo término de exploración. De hecho, será el valor del hiperparámetro 𝜉 el que determine la proporción de cada una de estas tendencias. Confidence Bound Criteria (UCB (Upper Confidence Bound) or LCB (Lower Confidence Bound)) En este caso se persiguen aquellos puntos que pudieran tener los valores más altos. 𝑈𝐶𝐵(𝒙)=𝜇(𝒙)+𝜅𝜎(𝒙),𝑐𝑜𝑛 𝜅≥0 (3-17) En el caso de escoger un valor 𝜅≤0, se estaría utilizando el LCB, equivalente a UCB para el caso de que se quiera minimizar la función objetivo. Caso Real. Resultados 26 26 La ejecución de cada experimento se divide en dos etapas. En la primera se lleva a cabo la configuración del mismo en base a los valores con los que se ha llamado a la función; de esta forma se definen la función de test, la función de adquisición, el kernel y el valor del parámetro de trade-off que se van a utilizar, así como las funciones de error y de precisión, el punto inicial y los límites de definición de la función de test. En la segunda etapa se lleva a cabo la implementación de la Optimización Bayesiana tal y como se indica en la Figura 4-4. En esta segunda etapa destacan dos aspectos importantes. El primero de ellos hace referencia a la estandarización de los datos obtenidos durante el muestreo. Este procedimiento es un requerimiento común de muchos algoritmos de machine learning, que basan su funcionamiento en que los datos sigan una distribución normal estándar, es decir, una distribución gaussiana con media cero y varianza unitaria. Esto no es más que realizar la siguiente operación sobre el valor de la muestra: 𝑦𝑠𝑡𝑎𝑛𝑑𝑎𝑟𝑑=𝑦𝑚𝑢𝑒𝑠𝑡𝑟𝑎−𝜇 𝑠 donde 𝜇 y 𝑠 son la media y la desviación estándar del conjunto de datos. Dado que el modelo se irá formando en base a los valores estandarizados, será necesario llevar a cabo la transformación inversa para calcular los indicadores de validez, de forma que dicho cálculo se haga con los valores reales equivalentes y no con los valores estandarizados. En segundo lugar, a la hora de definir el proceso gaussiano, se ha incluido un parámetro denominado alpha que modela el nivel de ruido en las muestras; matemáticamente, se trata de un valor constante añadido a la diagonal del kernel. Como es lógico, a mayor nivel de ruido, peor serán los resultados obtenidos. Sin embargo, ha sido necesario utilizar dicho parámetro por una cuestión de robustez computacional. En la expression (3-7) se vio que para actualizar los valores de media y varianza del proceso gaussiano es necesario llevar a cabo la inversión de la matriz K, proceso que puede dar lugar a inestabilidad numérica. Esta inestabilidad resulta mayor cuanto mayor sea el número de condición de la matriz, definido como la relación entre el mayor y el menor de sus autovalores. Teniendo en cuenta que los valores extremos de dichos autovalores aparecen cuando los valores de la función están muy correlacionados, resulta apropiado añadir el citado parámetro para prevenir problemas de inestabilidad numérica. Esto da lugar a un nuevo trade-off entre la robustez computacional y la precisión del modelo obtenido. Tanto el Código implementado en Python como los resultados presentados posteriormente pueden encontrarse en el repositorio de GitHub https://github.com/pjdg/TFG-Opt.-Bayesiana. 27 4.3 Resultados En esta sección se exponen los resultados de los experimentos realizados bajo las condiciones descritas anteriormente. Empleando las tres funciones de test anteriores y los kernels y funciones de adquisición que se presentaron en la sección 3 se han definido un total de 18 combinaciones posibles para las que obtener los indicadores de validez anteriores. Adicionalmente, cada una de estas configuraciones se ha ejecutado para varios valores del parámetro que controla el trade-off exploración/explotación en cada una de las funciones de adquisición. Más concretamente, el valor del parámetro de trade-off toma los valores 0.01, 0.1 y 0.2. Por otra parte, para los valores de los hiperparámetros que definen los kernels empleados se ha escogido 𝜈=1.5 y un valor de l igual al 10% de rango del mapa sobre el que se estudia la test function. Cabe destacar que en el caso de 𝑙 este únicamente será su valor inicial, ya que su valor será optimizado en cada iteración para mejorar el ajuste del kernel. Por otro lado, con el objetivo de obtener una estadística de los resultados obtenidos, se han realizado 30 réplicas de cada simulación. Esto se debe a que el algoritmo de optimización utilizado para optimizar la función de adquisición y obtener el siguiente punto a muestrear tiene una componente aleatoria, por lo que el conjunto de puntos muesteados varía de una ejecución a otra, y con él el valor de los indicadores finales. Para cada una de las diferentes configuraciones se ha construido una tabla resumen con los valores de error y desviación evaluados cada 10 iteraciones hasta un total de 50, siendo estos valores la media de los datos correspondientes a cada una de las réplicas simuladas. Además, la distribución del conjunto de estos datos podrá verse en las gráficas de caja y bigotes que acompañan a cada una de las tablas. También se expone un conjunto de gráficas que muestran la evolución de los mapas de error y de desviación conforme aumenta el número de iteraciones, así como los puntos que han sido muestreados sobre dichos mapas. Carece de sentido presenter el evolutivo para cada una de las réplicas simuladas, por lo que para cada experimento se muestra un único conjunto de mapas que ha sido escogido bajo criterio de mejor ajuste a la distribución de datos de dicho experimento. Caso Real. Resultados 28 28 Experimento 1: Branin - RBF - Probability of Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 44,07955 146,2179 20 27,76165 91,21267 30 22,93526 58,90926 40 18,78967 45,89285 50 16,07038 34,56775 0.1 10 34,55251 109,9385 20 19,92341 45,85731 30 10,27441 6,427224 40 8,155726 4,993284 50 7,134203 4,230362 0.2 10 37,88151 74,97493 20 16,95405 20,14319 30 9,648402 5,946918 40 8,156001 4,543914 50 6,351028 4,525501 Tabla 4-1. Resumen de resultados para el experimento 1 Figura 4-5. Resultados obtenidos para el experimento 1 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 29 Figura 4-6. Evolución de los mapas de error y precisión durante el experimento 1 Caso Real. Resultados 30 30 Experimento 2: Branin - RBF – Expected Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 29,95152 76,61173 20 16,00498 25,29415 30 9,953684 5,367825 40 8,124547 4,639944 50 6,540404 4,396046 0.1 10 31,40159 80,24767 20 13,99115 15,37785 30 9,207279 5,182346 40 7,478193 4,394321 50 6,373788 4,342641 0.2 10 30,8151 57,56132 20 13,33601 10,7587 30 8,783677 5,416258 40 7,4899 4,286748 50 6,32555 4,329187 Tabla 4-2. Resumen de resultados para el experimento 2 Figura 4-7. Resultados obtenidos para el experimento 2 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 31 Figura 4-8. Evolución de los mapas de error y precisión durante el experimento 2 Caso Real. Resultados 32 32 Experimento 3: Branin - RBF - Upper Confidence Bound 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 45,36464 94,01446 20 50,63141 122,5875 30 55,3247 139,4896 40 55,07298 143,3605 50 55,07749 143,4042 0.1 10 47,68209 121,3829 20 60,73948 155,552 30 60,766 155,843 40 60,76506 155,7414 50 60,71696 155,7442 0.2 10 53,97435 132,7817 20 60,95523 148,2527 30 60,75523 150,3981 40 60,76647 150,5017 50 60,61323 150,2205 Tabla 4-3. Resumen de resultados para el experimento 3 Figura 4-9. Resultados obtenidos para el experimento 3 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 33 Figura 4-10. Evolución de los mapas de error y precisión durante el experimento 3 Caso Real. Resultados 34 34 Experimento 4: Branin - Màtern kernel - Probability of Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 65,68003 114,9414 20 48,1067 88,85316 30 32,06024 49,85593 40 25,84069 38,19654 50 23,22447 31,03117 0.1 10 39,35278 95,39939 20 16,79825 33,95207 30 8,929918 13,91978 40 6,330839 9,551916 50 4,740604 7,698958 0.2 10 45,72681 69,4166 20 12,82491 25,8767 30 8,440863 11,97793 40 5,880425 8,977626 50 4,377671 7,370899 Tabla 4-4. Resumen de resultados para el experimento 4 Figura 4-11. Resultados obtenidos para el experimento 4 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 35 Figura 4-12. Evolución de los mapas de error y precisión durante el experimento 4 Caso Real. Resultados 42 42 Experimento 8: Schwefel - RBF – Expected Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 277,9958 257,0411 20 262,367 260,6042 30 255,068 245,71 40 249,168 238,0689 50 243,0912 229,3364 0.1 10 278,0818 263,9036 20 261,3414 259,8486 30 252,9058 247,1717 40 245,5073 235,312 50 236,5413 226,061 0.2 10 276,8218 264,2663 20 259,7398 265,7426 30 252,9667 253,0179 40 243,573 238,8563 50 236,4514 225,0281 Tabla 4-8. Resumen de resultados para el experimento 8 Figura 4-19. Resultados obtenidos para el experimento 8 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 43 Figura 4-20. Evolución de los mapas de error y precisión durante el experimento 8 Caso Real. Resultados 44 44 Experimento 9: Schwefel - RBF - Upper Confidence Bound 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 286,8533 262,0136 20 287,0928 279,7611 30 285,0387 278,9623 40 285,6591 286,4022 50 284,0568 287,5619 0.1 10 276,9246 272,9275 20 274,8248 290,7182 30 274,703 287,4644 40 274,7673 284,4386 50 274,7472 282,7443 0.2 10 270,3407 281,9114 20 268,2504 288,2165 30 268,189 288,8565 40 267,4405 288,5387 50 266,5523 292,4716 Tabla 4-9. Resumen de resultados para el experimento 9 Figura 4-21. Resultados obtenidos para el experimento 9 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 45 Figura 4-22. Evolución de los mapas de error y precisión durante el experimento 9 Caso Real. Resultados 46 46 Experimento 10: Schwefel - Màtern kernel - Probability of Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 296,283 279,5525 20 291,5839 294,0639 30 288,4264 296,199 40 288,0077 293,9555 50 287,7507 289,9383 0.1 10 284,406 296,2679 20 273,4523 272,5003 30 270,2511 257,4923 40 267,7147 245,4775 50 268,9378 233,7598 0.2 10 275,6236 267,2042 20 270,4353 255,7418 30 268,5634 242,6891 40 264,9959 237,2975 50 263,6105 229,2345 Tabla 4-10. Resumen de resultados para el experimento 10 Figura 4-23. Resultados obtenidos para el experimento 10 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 47 Figura 4-24. Evolución de los mapas de error y precisión durante el experimento 10 Caso Real. Resultados 48 48 Experimento 11: Schwefel - Màtern kernel - Expected Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 274,161 264,5666 20 267,2623 262,5879 30 266,3991 247,5796 40 265,9628 235,6715 50 260,0571 226,1262 0.1 10 276,6291 262,1144 20 265,9611 261,6339 30 261,9644 251,4827 40 257,8072 241,1697 50 254,3258 230,2749 0.2 10 276,532 268,8768 20 269,6792 257,957 30 265,4016 249,3654 40 259,3011 237,1734 50 251,9166 230,4637 Tabla 4-11. Resumen de resultados para el experimento 11 Figura 4-25. Resultados obtenidos para el experimento 11 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 49 Figura 4-26. Evolución de los mapas de error y precisión durante el experimento 11 Caso Real. Resultados 50 50 Experimento 12: Schwefel - Màtern kernel - Upper Confidence Bound 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 301,1965 211,8157 20 301,854 220,423 30 301,9836 236,0428 40 303,6296 240,4566 50 304,1494 240,683 0.1 10 303,6195 238,3356 20 305,3599 245,2556 30 305,518 243,0707 40 305,9632 242,4683 50 306,2622 242,415 0.2 10 294,5107 248,3997 20 295,6405 254,9291 30 296,7841 262,1055 40 296,9488 262,1932 50 297,0962 262,3435 Tabla 4-12. Resumen de resultados para el experimento 12 Figura 4-27. Resultados obtenidos para el experimento 12 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 51 Figura 4-28. Evolución de los mapas de error y precisión durante el experimento 12 Caso Real. Resultados 58 58 Experimento 16: Styblinski-Tang - Màtern kernel - Probability of Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 56,97751 148,7337 20 47,73169 139,7696 30 45,25832 131,4879 40 44,73806 126,6398 50 44,01831 124,1064 0.1 10 61,90941 134,459 20 60,58832 95,04539 30 35,30998 46,06423 40 12,03501 20,92078 50 7,438777 14,97738 0.2 10 87,03103 156,5593 20 45,63212 99,26978 30 18,66724 31,37681 40 9,802387 17,47103 50 6,119003 13,07478 Tabla 4-16. Resumen de resultados para el experimento 16 Figura 4-35. Resultados obtenidos para el experimento 16 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 59 Figura 4-36. Evolución de los mapas de error y precisión durante el experimento 16 Caso Real. Resultados 60 60 Experimento 17: Styblinski-Tang - Màtern kernel - Expected Improvement 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 63,97136 105,2279 20 50,32912 80,36717 30 19,98814 29,80382 40 9,719155 18,17172 50 5,669146 13,45697 0.1 10 60,43697 106,3266 20 49,20846 79,8255 30 17,36313 27,2876 40 9,11913 17,22617 50 5,672913 12,92843 0.2 10 60,24561 108,7522 20 48,9526 80,16455 30 16,42493 26,31698 40 9,058595 17,01159 50 5,512996 12,74784 Tabla 4-17. Resumen de resultados para el experimento 17 Figura 4-37. Resultados obtenidos para el experimento 17 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 61 Figura 4-38. Evolución de los mapas de error y precisión durante el experimento 17 Caso Real. Resultados 62 62 Experimento 18: Styblinski-Tang - Màtern kernel - Upper Confidence Bound 𝜉 n 𝑒𝑐𝑚 𝜎𝑐𝑚 0.01 10 46,3315 107,2135 20 105,2183 23,47403 30 75,16282 63,87796 40 75,30469 63,42633 50 75,35785 63,32431 0.1 10 130,0195 67,23454 20 130,0185 67,03853 30 130,0227 66,96965 40 130,0276 66,9196 50 130,0319 66,88169 0.2 10 80,25931 86,54647 20 80,0187 86,53393 30 79,81513 86,53393 40 79,81609 86,55375 50 79,8162 86,56698 Tabla 4-18. Resumen de resultados para el experimento 18 Figura 4-39. Resultados obtenidos para el experimento 18 a las n iteraciones. De izquierda a derecha: 𝜉=0.01, 𝜉=0.1, 𝜉=0.2 63 Figura 4-40. Evolución de los mapas de error y precisión durante el experimento 18 Caso Real. Resultados 64 64 4.4 Discusión de los resultados Con el objetivo de obtener conclusiones adecuadas respecto a las diferentes configuraciones para las que se han presentado los resultados anteriores, se realizará un análisis de estos en base a los diferentes elementos configurables que se han empleado -test function, kernel, función de adquisición y parámetro de trade-off- a modo de descenso desde el elemento más general al más concreto. En este recorrido será de gran ayuda la figura 4-41, que reúne los resultados genéricos de cada experimento a modo de resumen. En primer lugar, los modelos obtenidos cuando se intenta aproximar la función de Schwefel no resultan acertados en ninguno de los casos propuestos: en todos ellos el error alcanza valores entre 200 y 300, lo que hace que los modelos resulten inadmisibles para la finalidad perseguida. Aún más, el incremento en el número de iteraciones no conlleva una mejora sustancial en el modelo (llegando incluso a empeorar el error en el caso del experimento 11) con una disminución del error de un 15% (medido entre 10 y 50 muestras) en el caso más favorable. Además, a medida que se toman más muestras también se aprecia un aumento en la incertidumbre del modelo, por lo que el conjunto de las dos tendencias descritas lleva a pensar que la calidad del modelo empeorará a medida que aumente el número de iteraciones. Sin embargo, el panorama cambia cuando se pretenden aproximar las dos funciones de test restantes. En ambos casos se llega a obtener buenos resultados con todas las configuraciones: salvando aquellos casos en los que se emplea UCB como función de adquisición, en los que se incidirá más adelante, el error final obtenido más alto ronda un valor de 6, mientras que en los casos más favorables baja a 4 para la función de Branin y a 2 para la función de Styblinski-Tang. Además, en estos casos la evolución a medida que aumenta el número de muestras tomadas resulta mucho más favorable tanto para el error como para la incertidumbre del modelo, logrando disminuciones de error del 91% y del 96% en los casos para Branin y para Styblinski-Tang, respectivamente; en definitiva, el aumento en el número de iteraciones hace que que el modelo se vuelva más exacto y preciso. Figura 4-41. Resumen de los resultados de error para n iteraciones. 65 La polarización de resultados entre la pareja formada por estas dos últimas funciones y la función de Schwefel se debe fundamentalmente a la elección del kernel; como ya se explicó en la sección 3, las consecuencias de esta decisión predominan en el resultado del algoritmo. Los kernels Mátern y RBF resultan ser más suaves de lo que lo es la función de Schwefel, lo que dificulta enormemente que el modelo pueda adaptarse a una función de este tipo. Además, el hecho de que dicha función presente tal número de máximos y mínimos (Figura 4-2) dificulta enormemente que un algoritmo cuyo próposito original es la optimización lleve a cabo un proceso de exploración más exahustivo, tendiendo a tomar muestras en zonas más concretas; esto impide la recopilación de información que le permita mejorar la predicción de la función. Por su parte, las funciones de Branin y de Styblinski-Tang resultan ser mucho más suaves y semejantes entre sí, lo que resulta coherente con la mejora y la similitud presentes en los resultados obtenidos en estos casos. Cabe destacar que todos estos resultados no deben llevar a la conclusión de que la Optimización Bayesiana no es apta para el modelado de funciones como la de Schwefel. Podrían conseguirse buenos resultados en tanto que la elección del kernel empleado resulte más acorde con esta función. Sin embargo, en este trabajo no se ha querido seguir este camino dado que hacerlo entraría en conflicto con la realidad subyacente del caso estudiado: en la monitorización del lago Ypacaraí el modelo del lago se considera completaente desconocido, por lo que resulta de interés medir qué tan buenos resultados pueden obtenerse enfrentando un kernel predeterminado a diferentes entornos, y no a la inversa. Sin ánimo de incidir más sobre la importancia del kernel en el desempeño del algoritmo, puede apreciarse que en el caso de la función de Schwefel esta elección no juega un papel muy relevante, mientras que los dos casos restantes si presentan ciertas diferencias. Si bien en el caso de la función de Branin se obtienen mejores resultados cuando se emplea el kernel Màtern (con un valor del error de 4 frente a 6 si se emplea el RBF), al intentar modelar la función de Styblinski-Tang es preferible usar el kernel RBF (obteniendo un error de 2 frente al 6 del Màtern). En todos los demás casos se comprueba que el uso del RBF resulta favorable frente al uso del Màtern. Esta pequeña contraposición de resultados entre los experimentos con las funciones de Styblinski-Tang y de Branin se debe, una vez más, al grado de adecuación del kernel a la función de test. En el primer caso, la función presenta un comportamiento más similar al exponencial, principalmente en zonas cercanas a la frontera definida, por lo que es lógico que el RBF obtenga mejores resultados. Sin embargo, la función de Branin se aleja ligeramente de este comportamiento, lo que hace que este kernel empeore su desempeño siendo superado por el Màtern (más flexible). Independientemente del kernel utilizado, la gráfica resumen anterior permite comprobar que, en líneas generales, los resultados obtenidos al emplear las funciones de adquisición PI y EI no presentan una diferencia sustancial entre ellos, mientras que sí difieren con respecto a los datos obtenidos al utilizar UCB, siendo estos resultados significativamente peores y apenas variantes con el número de muestras tomadas. Esto se debe a que la función de adquisición UCB tiene un enfoque de explotación pura, haciendo que el muestreo se realice de manera más focalizada (como puede comprobarse en la evolución de los mapas para estos casos), con la consecuente falta de información debido a una exploración deficiente. Caso Real. Resultados 66 66 En el caso de las funciones de Branin y de Styblinski-Tang, que resultan más favorables para la exploración, el error final alcanzado aumenta un orden de magnitud cuando se pasa de usar PI o EI a usar UCB como funciones de adquisición. El caso de la función de Schwefel resulta algo más particular, ya que su abundancia de máximos y mínimos difculta enormemente la exploración. Si bien puede comprobarse en las gráficas presentadas anteriormente cómo se mantiene la tendencia de mayor exploración al usar PI o EI, las muestras suelen tomarse en los puntos donde la función toma valores extremos, por lo que el error cometido por el modelo no se reduce notablemente con respecto al alcanzado al usar UCB. Por tanto, queda demostrada la importancia de una correcta exploración para conseguir buenos resultados; de ahí la relevancia del balance entre la tendencia innata del algoritmo hacia la explotación y la búsqueda necesaria de la exploración, en el que el parámetro de trade-off 𝜉 toma relevancia. La Figura 4-40 se ha construido tomando los resultados más favorables de cada caso, omitiendo la información relativa a la variación del parámetro 𝜉 pero que sí está presente en las tablas de resultados expuestas anteriormente. Por una parte, puede comprobarse que, aunque se ha dicho que las funciones de adquisición PI y EI no presentan una diferencia sustancial en cuanto a los resultados obtenidos, sí que lo hacen para el valor más bajo del parámetro de trade-off, es decir, cuando este apenas tiene efecto sobre las propias funciones. En estos casos se aprecia la mejora proporcionada por la función EI frente a PI, consiguiendo una diminución del error del final de hasta el 78% en el caso de la función de Branin y de un 95% para la función de Styblinski-Tang. Por otro lado, el efecto de este parámetro sobre la función UCB puede apreciarse en los mapas de evolución, teniendo lugar una dispersión de puntos de muestreo mayor a medida que el valor del parámetro aumenta: valores mayores favorecen la exploración, contrarrestando en cierta medida la fuerte tendencia de explotación de la función UCB. 67 5 IMPLEMENTACIÓN REAL Hasta el momento se ha llevado a cabo un estudio y análisis de la Optimización Bayesiana como técnica de modelado de un entorno desconocido. Es importante conocer el método y su alcance, pero también lo es conocer su implementación en el caso práctico real, por lo que resulta importante llevar a cabo una breve descripción de cuál es concretamente el papel que juega la Optimización Bayesiana en el sistema del ASV [28]. El sistema que gobierna el movimiento del vehículo se trata de un GNC (Guidance, Navigation and Control), cuyo esquema queda representado en la Figura 5-1. Dentro de este conjunto, la Optimización Bayesiana pertenece al bloque de orientación, encargado de decidir el punto de destino y la ruta a seguir hasta este. A grandes rasgos, el sistema presenta un comportamiento cíclico en el que el sistema de navegación y los sensores aportarán al sistema de orientación la información necesaria para decidir el movimiento que será ejecutado por el sistema de control del vehículo. El sistema de posicionamiento deberá tener una precisión aproximadamente igual a la resolución del mapa que se haya definido para el entorno, de manera que en todo momento el vehículo quede localizado en uno de los puntos que conforman dicho mapa. Por otra parte, el sistema de orientación debe incluir un algoritmo de planificación de trayectoria local que que lleve a cabo la evasión de obstáculos mientras el vehículo recorre el camino descrito por el algoritmo de planificación global, en este caso la Optimización Bayesiana. Figura 5-1. Esquema del sistema GNC del ASV. Fuente: [28]