**IAIC 2025-26** Optimización  !!!side:1 Veremos en la segunda parte del curso que esto es especialmente cierto en el campo del Aprendizaje Automático (Machine Learning), donde un porcentaje muy amplio de los problemas se ven bajo el cristal de la Optimización. La optimización de funciones es un tema de gran importancia, tanto teórica como práctica, en el campo de las matemáticas en general (con aplicaciones en todas las ingeinerías) y de la IA en particular. Los problemas de optimización aparecen en muchos contextos diferentes, bajo el escenario común de optimizar una medida objetivo con respecto a un conjunto de parámetros $^1$. Podemos considerar distintas clasificaciones dentro de la optimización, según atendamos a unas características u otras: * En cuanto al tipo de variables que usa la función objetivo, podemos diferenciar entre la optimización **discreta** y la **continua**. En la optimización discreta, algunas o todas las variables de la función objetivo son discretas, es decir, sólo asumen valores de un conjunto finito de valores o, como mucho, un conjunto disperso de valores similar a los números naturales. Dos campos importantes de la optimización discreta son la **Optimización Combinatoria** y la **Programación Entera**. La optimización continua suele aprovecharse de propiedades adicionales, y no solo de la continuidad, como es la diferenciabilidad, y podemos encontrar métodos que hacen unos de las primeras derivadas para dar algoritmos de existencia y aproximación de los óptimos locales (que es un primer paso para encontrar los óptimos globales). * Además, podemos atender a una clasificación en función de la presencia de restricciones sobre las variables independientes de la función objetivo (de forma similar a como vimos en CSP). Las restricciones pueden ser **duras**, que deben satisfacer obligatoriamente las variables independientes, o **blandas**, en las que se penaliza la función objetivo si no se satisfacen las restricciones (el valor de la penalización suele ir en función de la medida en que no se satisfagan), pero se admiten soluciones que no las verifiquen. !!!side:2 Si, además de ser una expresión analítica, es diferenciable respecto a los parámetros, entonces hay toda una batería de metodos de optimización basados en el gradiente de la expresión. * En la elección del método de optimización también influye el coste computacional de la evaluación de la función objetivo. Si la función objetivo se da como una expresión analítica, normalmente es posible realizar muchas evaluaciones de la función, lo que hace que el problema de optimización sea más manejable $^2$. Pero si la función objetivo es costosa desde el punto de vista computacional (como el resultado de un experimento que hay que hacer en un laboratorio de química o de una simulación lenta), la elección del método de optimización es más limitada y también más importante. En este tema vamos a comenzar viendo la relación existente entre **Búsquedas** y **Optimización**, con el fin de usar los algoritmos que desarrollemos en un campo para resolver los problemas que se pueden plantear desde el otro. Como ya vimos en el tema de búsquedas, podemos hacer uso del proceso de optimización de heurísticas como estrategias para buscar en el dominio de forma *inteligente*. !!!note Aunque en el tema anterior hemos indicado que es normal añadir la restricción de que el Espacio de Estados sea finito, algo conveniente para los métodos seguidos en las búsquedas clásicas, ahora no es estrictamente necesario, y los métodos que veremos aquí también se pueden aplicar en casos en los que el Espacio de Estados es infinito, incluso no numerable (como los espacios euclídeos, $\mathbb{R}^n$). Entre los métodos heurísticos más habituales para optimización podemos encontrar la **evolución diferencial**, la **computación evolutiva** (incluyendo, por ejemplo, los **algoritmos genéticos** y la **programación genética**), el **templado simulado**, y la **optimización social** (por ejemplo, basada en **enjambres de partículas** y en **colonias de hormigas**). En este tema nos centraremos únicamente en un par de métodos con el fin de complementar las aproximaciones vistas en los temas anteriores: el **Templado Simulado** y la **Optimización por Enjambre de Partículas** (**PSO**). Dejamos fuera el famoso método de **Algoritmos Genéticos**, no porque carezca de interés, sino porque es probable que el alumno ya lo haya estudiado en otras asignaturas o lo vaya a hacer en el futuro. La cantidad de variantes existentes plantea la cuestión fundamental de qué método utilizar cuando se presenta un problema de optimización. Por desgracia, no hay una respuesta general a esta pregunta, ya que **no existe el mejor algoritmo de optimización**. Una realidad que viene formalizada en los **Teoremas No-Free-Lunch**: !!!note: Teoremas No-Free-Lunch Los **Teoremas No-free-lunch** son resultados en optimización y complejidad computacional que implican que, para ciertos tipos de problemas, el coste computacional de encontrar una solución cuando se promedia sobre todos los problemas de una clase es el mismo para cualquier método de solución. En otras palabras, **no existe ningún método de optimización que supere a todos los demás métodos en todos los problemas**. Cualquier mejora del rendimiento en una clase de problemas siempre se ve compensada por una disminución del rendimiento en otra clase de problemas. Los métodos de optimización especializados tienen el mejor rendimiento para resolver una determinada clase de problemas, pero suelen dar resultados muy pobres cuando se aplican a clases distintas. En cambio, aunque puede haber algoritmos que den buenos resultados en diversas clases de problemas, siguen siendo superados por los diversos métodos especializados que se han diseñado para cada una de ellas. # Búsquedas por Optimización Hasta ahora hemos visto métodos de resolución de problemas en los que el objetivo ha sido buscar un camino en un espacio de estados, ya sea porque estamos interesados en alcanzar un estado final (búsqueda de estados simple), o porque es el propio camino el que nos interesa (planificación). Por desgracia, los métodos vistos para dar estas soluciones no son siempre factibles, por dos razones principalmente: * Porque este planteamiento resulte **demasiado artificial** y lejano al problema, por ejemplo, porque no sea natural definir operadores de transición, o porque no sea una tarea sencilla definir una función heurística adecuada para la búsqueda (que sea consistente y admisible, por ejemplo, cuando usamos $A^*$). * Porque no hay pocas garantías de hallar una solución óptima debido a que el **tamaño del espacio de búsqueda es demasiado grande**, y en este caso nos podríamos conformar con una solución que podamos considerar **suficientemente** buena. Si estamos ante un caso así, podemos optar por dos vías: 1. Si es relativamente asequible encontrar una solución inicial, aunque no sea demasiado buena, podemos cambiar nuestra tarea de buscar un camino entre estados a la de **encontrar una mejora sucesiva** de soluciones no óptimas en el espacio de soluciones. 2. Si no es posible encontrar una solución inicial (aunque no sea suficientemente buena) podemos ver si, partiendo de un estado no válido, podemos ir construyendo una sucesión de estados que sí se aproximen a una solución. En este cambio de perspectiva, pues, lo que suponemos es que podemos navegar dentro del espacio de estados, realizando operaciones sobre ellos que pueden ayudar a mejorarlos. Estas operaciones no tienen por qué tener una representación en el mundo real, como si fueran las reglas de un juego conocido a priori, sino que son simplemente formas de manipular las posibles soluciones para encontrar otras opciones. El coste de estas operaciones no será tan relevante, ya que lo que estamos buscando es la calidad de las soluciones, y no la forma de llegar a ellas (por tanto, los métodos que veamos aquí son más adecuados a búsquedas que a planificación). Al igual que en temas anteriores disponíamos de una heurística que informaba de la bondad de hacer una elección frente a otra, ahora necesitamos una función que, según algún criterio (dependiente del dominio), mida la calidad de las posibles soluciones y permita ordenarlas respecto a su bondad. La estructura general que tendremos en este tipo de aproximaciones es similar a la que ya vimos en las búsquedas en estados de espacios: 1. Una **propuesta de solución inicial**. 2. **Operadores de transformación de la solución**, que permiten ir generando o modificando las propuestas de solución actuales. No tienen coste explícito asociado, pero se suelen buscar muy ligeras computacionalmente para poder hacer muchas pruebas si fuera necesario. 3. Una **función de calidad** (a veces llamada **función de fitness**), que mide lo buena que es una propuesta de solución y permite dirigir la búsqueda. A diferencia de las funciones heurísticas anteriores, en este caso esta función no indica cuánto falta para encontrar la solución que buscamos (y tampoco deben verificar condiciones adicionales como admisibilidad o consistencia), sino que proporciona una forma de comparar la calidad entre soluciones para indicarnos qué transformaciones pueden ser las más adecuadas. Cuando se aplican estos algoritmos de aproximación a funciones de calidad que pueden tomar valores continuos, saber cuándo hemos de detener la búsqueda dependerá totalmente del problema y del algoritmo, que tendrá que decidir cuándo ya no merezca la pena seguir avanzando (por ejemplo, porque las mejoras no compensan el esfuerzo computacional adicional). !!!note: Optimización vs. Búsqueda Lo que hace una función de fitness es convertir el problema de búsqueda en un problema de optimización. En muchos casos interesantes, como los que a menudo encontramos en problemas de Inteligencia Artificial, no disponemos de funciones de calidad con una expresión analítica conocida o que verifiquen las condiciones necesarias para aplicar algoritmos analíticos de búsqueda de óptimos en funciones, por lo que tenemos que explorar la función de calidad evaluándola en los vecinos de la solución actual, es decir, haciendo uso únicamente de información local durante el proceso. !!!warn:Minimización o Maximización En lo que sigue, y guiados por el uso de la optimización para abordar los problemas de búsqueda con heurísticas (donde la búsqueda se traduce en minimizar la función heurística), vamos a estudiar el caso de la **minimización** de funciones. En todo caso, transformar una minimización en una maximización (o viceversa) es algo directo, ya sea adaptando trivialmente los algoritmos o teniendo en cuenta que: $$\max \{f(x): x\in D\} = \min \{-f(x):\ x\in D\}$$ # Templado Simulado Muchas veces los algoritmos de optimización se quedan atascados enseguida en óptimos locales, lejos de los valores de un óptimo global. Para evitar esta convergencia temprana se han diseñado algoritmos que tienen como tarea escapar de esas soluciones subóptimas con el fin de poder explorar una maor sección del espacio y mejorar las soluciones ofrecidas. El algoritmo del **Templado Simulado** (también conocido como **Simulated Annealing**, o **Recocido Simulado**) es uno de ellos y está inspirado en un fenómeno físico que se observa en el templado de metales y en la cristalización de disoluciones: !!!note: Templado Simulado en el Mundo Físico Todo conjunto de átomos o moléculas tiene un estado de energía que depende de la temperatura del sistema. A medida que el sistema se enfría, va perdiendo energía hasta que se estabiliza debido a los enlaces que se van formando entre los elementos del sistema (y que, cuando la temperatura es adecuadamente alta, no tienen fuerza suficiente para mantenerse y se crean y destruyen de forma continua). El fenómeno físico que se observa es que, dependiendo de cómo se realiza el enfriamiento, el estado alcanzado de energía final es muy diferente. Por ejemplo, si una barra de metal se enfría demasiado rápido la estructura cristalina a la que se llega al final está poco cohesionada (un número bajo de enlaces entre átomos, o mal distribuidos) y, consecuentemente, se rompe con facilidad. Si el enfriamiento se realiza más lentamente, el resultado final es totalmente diferente, el número de enlaces entre los átomos es mayor y está mejor distribuido, lo que proporciona una barra mucho más difícil de romper. Durante este enfriamiento controlado se observa que la energía del conjunto de átomos no disminuye de manera constante, sino que a veces la energía total puede ser mayor que en un momento inmediatamente anterior (y, si es necesario, se aporta energía artificalmente, como hace el herrero en el templado del hierro por medio de la fragua). Esta circunstancia hace que el estado final que se alcanza bajo estas condiciones controladas sea mejor que cuando el enfriamiento se realiza rápidamente y de forma natural. Inspirados en este fenómeno físico podemos idear un algoritmo que permite, para ciertos problemas, obtener soluciones mejores que las proporcionadas por los algoritmos de búsqueda vistos hasta ahora. En este algoritmo, el estado siguiente a explorar no será siempre el mejor vecino (que es lo que haría un algoritmo voraz), sino que se elige con cierta aleatoriedad de entre todos los vecinos (mejores o peores). Una ventaja de este algoritmo es que no hay por qué generar todos los sucesores de un estado (algo que puede ser muy costoso, o imposible en el caso continuo), basta elegir un sucesor al azar y decidir si continuamos por él o no (por lo que es un algoritmo estocástico). Además, solo requiere dos operaciones sencillas: la elección de un punto de partida aleatorio, y una operación que calcula un único un posible sucesor, y puede aplicarse a problemas de optimización y búsqueda tanto continuos como discretos. ![](./img/simann.jpg align="right" width=50%)La nomenclatura del algoritmo de Templado Simulado transporta la analogía física al problema que se quiere solucionar: * Se denomina **función de energía** a la función que mide la calidad de una solución (recordemos que en el problema físico se pretende minimizar la energía del sistema). * Un parámetro, denominado **temperatura**, controla el funcionamiento del algoritmo. * La temperatura y la diferencia de calidad (incremento de energía) entre el estado actual y el sucesor generado determinan si se acepta la nueva propuesta o no: * Si el sucesor es mejor, se acepta. * Cuanta mayor sea la temperatura, más probabilidad habrá de que un sucesor peor sea aceptado. * La elección de un estado peor dependerá de su diferencia de calidad con el estado actual, cuanta más diferencia, menos probabilidad de aceptarlo. * Si el sucesor generado es peor y no se acepta, se genera otro sucesor. !!!side:3 Aunque hay algunos resultados teóricos acerca de la influencia de estos parámetros, en general hay que decidir experimentalmente cuál es la temperatura inicial más adecuada y también la forma más adecuada de hacer que vaya disminuyendo. * Se sigue una estrategia de enfriamiento: cada cierto número de iteraciones, el sistema se enfría un cierto ratio. El proceso para cuando la temperatura llega a un umbral mínimo prefijado. La elección de estos dos parámetros determina el comportamiento completo del algoritmo, y puede modificar su eficiencia de forma considerable $^3$. ## Algoritmo La representación en forma de pseudocódigo de esta idea es la siguiente (no es la única posible representación, ya que hemos fijado varios de los grados de libertad que hemos comentado): ~~~~~~~~ Templado Simulado (T₀, s₀, N, γ < 1, ϵ > 0) T = T₀ while T > ϵ Repeat N: s₁ = Genera_vecino(s₀) ∆E = E(s₀) − E(s₁) if ∆E > 0 s₀ = s₁ else con probabilidad exp(∆E/T): s₀ = s₁ T = T ⋅ γ return s₀ ~~~~~~~~ Como en la analogía física: 1. Si el enfriamiento es muy rápido ($\gamma<<$), la temperatura bajará muy rápidamente, y el camino explorado será relativamente aleatorio, porque no tiene tiempo de comparar muchas opciones. 2. En caso contrario, la bajada de temperatura será más suave y la búsqueda será mejor. 3. Si $N$ es muy grande, localmente podrá tener más opciones de elegir vecinos en los que se va produciendo una mejora. El algoritmo favorece la disminución de la energía, $E$, ya que: 1. Si $\Delta E>0$ ($E(s_0)>E(s_1)$), entonces $s_1$ se convierte en el estado actual. 2. Cuando $\Delta E < 0$, se tiene que $e^{\Delta E/T} < 1$ y, de hecho: $$\displaystyle \lim_{T\to 0} e^{\Delta E/T}=0$$ por lo que el hecho de admitir nuevos estados que empeoran la búsqueda se hace cada vez más improbable a medida que avanza el algoritmo y la temperatura decae. Además, los estados que empeoran mucho la energía (es decir, aquellos que hacen que $\Delta E$ sea más negativo) se hacen menos probables que los que la empeoran menos. La función de energía tiene muchas similitudes con las heurísticas de algoritmos como $A^∗$, ya que de alguna forma mide lo buena, o cerca de una solución ideal, que está el estado actual, pero tiene una ventaja computacional, y es que no imponemos ninguna propiedad adicional como ser admisible, aunque sí se le supone un comportamiento bastante básico: 1. $E(s)\geq 0$, para cualquier estado, $s$, de nuestro espacio. 2. $E(s)$ mínimo si y solo si $s$ es un óptimo. 3. $E(s_1) < E(s_2)$, si $s_1$ es *mejor solución* que $s_2$. El Templado Simulado es especialmente bueno para problemas con un espacio de estados grande en los que el óptimo global está rodeado de muchos óptimos locales. Se puede utilizar también para problemas en los que encontrar una heurística discriminante es difícil (donde una elección aleatoria es tan buena como otra cualquiera), ya que en estos casos muchos otros algoritmos no tienen mucha información con la que trabajar y se quedan atascados enseguida, mientras que el templado simulado puede explorar un espacio de soluciones más grande y tiene mayor probabilidad de encontrar una solución buena. Su mayor problema es determinar los valores de los parámetros, y requiere una importante labor de experimentación que depende de cada problema. Los valores óptimos de estos parámetros varían con el dominio e incluso con el tamaño de la instancia concreta del problema. !!!code:Julia

El algoritmo anterior se puede implementar muy sencillamente en Julia sin necesidad de usar estructuras ni funciones auxiliares. En la implementación que veremos se hace uso de las siguientes variables y funciones: - `Energia::Function`: Función de energía/fitness a minimizar (toma estados). - `s0`: Estado inicial. - `N::Int`: Número de repeticiones por cada ciclo de temperatura, `T`. - `vecino::Function`: Función que devuelve un vecino del estado actual. - `T_min::Float64`: Valor mínimo de Temperatura. - `tasa_T::Float64`: Tasa de enfriamiento (medido como porcentaje de disminución de `T`). Es equivalente a $\gamma$ en el algoritmo original. - `ac_eq::Bool`: Indica si se acepta el cambio cuando $ΔE=0$. ~~~~ Julia linenumbers function SA(Energia::Function, s0, N::Int, vecino::Function, T_min::Float64, tasa_T::Float64, ac_eq::Bool) # Datos de partida T = 1.0 sbest = s0 Ebest = E0 = Energia(s0) # El algoritmo itera hasta alcanzar la T_min while T > T_min # Hacemos N pruebas para un T fijo for _ in 1:N s1 = vecino(s0) # Tomamos un vecino (generador) E1 = Energia(s1) # Calculamos su energía ΔE = E1 - E0 # y el incremento de energía # Si la energía disminuye, o es igual pero se acepta, o con prob∼T if ΔE ≤ 0 || (ac_eq && ΔE ==0) || rand() < T s0 = s1 # saltamos a ese estado E0 = E1 # Si el nuevo estado es el mejor if E0 < Ebest sbest = s0 # lo almacenamos Ebest = E0 end end end # Actualizamos T T = T * (1 - tasa_T / 100) end # Devolvemos el mejor estado encontrado, y su energía return sbest, Ebest end ~~~~ El algoritmo devuelve el mejor estado encontrado en el proceso de búsqueda, y el valor de `Energia` correspondiente a ese estado.

Veamos un ejemplo de aplicación a un caso ya estudiado anteriormente: el problema de las $N$ reinas. !!!code:Julia

Recordemos el enunciado del problema: en un tablero de ajedrez de tamaño $N\times N$, colocar $N$ reinas de manera que ninguna amenace a otra (recordamos que dos reinas se amenazan si están en la misma fila, columna, o diagonal). Una propuesta de representación sencilla consiste en un vector de $N$ posiciones donde la posición $i$ indica la fila de la reina $i$-ésima (es decir, por defecto, cada reina ocupa una columna distinta). Bajo esta representación, una propuesta de vecino es modificar al azar uno de los valores (también elegido al azar) del vector de $N$ posiciones: ~~~~ Julia linenumbers function vecino(s) N = length(s) s1 = copy(s) i = rand(1:N) # Elegimos una reina al azar s1[i]=rand(1:N) # y modificamos su posición (en la misma columna) return s1 end ~~~~ Podemos definir una función auxiliar que indica si hay una amenaza entre la reina que ocupa la columna $i$ y la que ocupa la $j$, con el fin de contar cómo de buena (cuántas amenazas hay) es una configuración: ~~~~ Julia linenumbers # ¿La reina i amenaza a la reina j en el estado s? function amenaza(i,j,s) eli = s[i] elj = s[j] return (eli == elj || abs(elj - eli) == abs(j - i)) ? 1 : 0 end ~~~~~ El número total de amenazas (que será nuestra función de energía), es simplemente contar el número total de amenazas en el estado. Téngase en cuenta que el mínimo de esta función de energía (cuando vale $0$) coincide exactamente con que el estado asociado sea solución del problema ($N$ reinas colocadas sin amenazas entre ellas): ~~~~ Julia linenumbers function amenazas(s) N = length(s) res = 0 for i in 1:N for j in (i+1):N res += amenaza(i,j,s) end end return res end ~~~~ A partir de estas definiciones podemos resolver ya el problema por medio del algoritmo de Templado Simulado. Para ello, vamos a partir de un estado inicial completamente aleatorio, y dejaremos que el algoritmo vaya mejorando su función de energía (`amenazas`) con el fin de aproximarlo a $0$. ~~~~Julia # Estado inicial aleatorio N = 8 s0 = map(x -> rand(1:N), 1:N) # Ejecutamos Templado simulado a partir del estado anterior solReinas, amReinas = SA(amenazas, s0, 1000, vecino, 0.0001, 1.0, true) ~~~~ Julia output ([3, 7, 2, 8, 6, 4, 1, 5], 0) ~~~~ Que representa la siguiente solución (también realizada en Julia, y que está en los ficheros que se distribuyen con el curso): 

Veamos ahora cómo resolver el Problema del Viajante de Comercio: !!!code:Julia

Dado un conjunto de $N$ ciudades conectadas entre sí (en una primera versión, podemos suponer que todas están conectadas con todas), encontrar una ruta que pase por todas y cada una de ellas una sola vez (lo que se denomina un *tour*), y con recorrido total mínimo. Vamos a suponer que las ciudades ocupan posiciones 2D dentro de $[0,1]\times[0,1]$. El primer paso es definir una estructura que permita representar un tour que visita todas las ciudades. Esta estructura consta de los siguientes campos: - `orden::Vector{Int}`: orden en el que se visita cada ciudad. - `D::Matrix`: Matriz de distancias. Se almacena para no tener que recalcularla cada vez (se podría haber puesto como una variable de acceso global, pero hubiera quedado más artifical el uso en el algoritmo). Y daremos algunas funciones auxiliares útiles: - `RandomRuta(D::Matrix)`: Crea un tour con un orden aleatorio de las ciudades. La matriz de distancias permanece igual.n que proporciona una ruta aleatoria inicial - `cost(tour::Tour)`: Calcula el coste total de un tour sumando las distancias entre ciudades sucesivas. Será usado como función de energía a minimizar. - `cambia_tour(tour::Tour)`: Propone un nuevo tour candidato a partir de otro invirtiendo el orden de las ciudades entre dos índices aleatorios. Será usado como función que calcula un estado (tour) vecino a partir de otro. ~~~~ Julia linenumbers struct Ruta orden::Vector{Int} D::Matrix{Float64} end # randcycle necesita el paquete Random using Random RandomRuta(D::Matrix) = Ruta(randcycle(size(D, 1)), D) function coste(ruta::Ruta) c = 0.0 for (k, c1) in enumerate(ruta.orden) c2 = ruta.orden[mod1(k + 1, length(ruta.orden))] c += ruta.D[c1, c2] end return c end function cambia_tour(ruta::Ruta) n = length(ruta.orden) orden = copy(ruta.orden) # El orden se invertirá entre los índices i1 e i2 (incluidos) i1 = rand(1:n) i2 = mod1(rand(1:n-1) + i1, n) # Debe asegurarse que i1 != i2 # y que i1 < i2 if i1 > i2 i1, i2 = i2, i1 end if i1 == 1 && i2 == n return Ruta(orden, ruta.D) end orden[i1:i2] = reverse(orden[i1:i2]) return Ruta(orden, ruta.D) end ~~~~ Para usar el Templado Simulado para resolver este problema, necesitamos dar un conjunto inicial de ciudades (aleatorias, por ejemplo), calcular la matriz de distancias entre ellas, y llamar a la función `SA` de forma adecuada: ~~~~ Julia # Conjunto de ciudades aleatorio. Coste mínimo desconocido num_ciudades = 50 ciudades = [transpose(rand(num_ciudades)); transpose(rand(num_ciudades))] ~~~~ Julia output 2×50 Matrix{Float64}: 0.443839 0.603455 0.802701 0.994229 … 0.829244 0.958505 0.465726 0.61283 0.277974 0.390222 0.194905 0.393513 0.12935 0.953327 ~~~~ Julia # Matriz de distancias (no cambia a lo largo del experimento) D = sqrt.((ciudades[1, :] .- ciudades[1, :]').^2 + (ciudades[2, :] .- ciudades[2, :]').^2) ~~~~ Julia output 50×50 Matrix{Float64}: 0.0 0.370952 0.422299 0.691079 … 0.443438 0.706141 0.3412 0.370952 0.0 0.228689 0.399506 0.253634 0.384903 0.689254 0.422299 0.228689 0.0 0.273554 0.0267467 0.303857 0.656232 0.691079 0.399506 0.273554 0.0 0.258196 0.0746572 0.924402 0.23829 0.132739 0.25967 0.494001 0.286158 0.493731 0.562476 0.619417 0.382285 0.198458 0.115846 … 0.176019 0.182819 0.826096 0.398544 0.309469 0.106633 0.344458 0.101215 0.389088 0.582322 0.449957 0.562282 0.750128 0.961782 0.776676 0.942858 0.698085 0.439365 0.621389 0.788725 1.01832 0.814632 1.00609 0.639433 0.306957 0.677784 0.704253 0.977778 0.721151 1.00126 0.178352 ⋮ ⋱ 0.55367 0.215975 0.201026 0.193968 0.206337 0.169109 0.834174 0.791049 0.447425 0.403772 0.173052 0.396495 0.10676 1.05854 0.546137 0.175873 0.311916 0.342741 0.327606 0.297348 0.86383 0.503739 0.590128 0.786843 0.988548 0.813534 0.965753 0.755771 0.124565 0.375662 0.354961 0.628438 … 0.371156 0.655359 0.317018 0.587822 0.342506 0.165541 0.121023 0.145365 0.176579 0.80584 0.509068 0.31344 0.0959972 0.211788 0.0699156 0.261448 0.715279 0.373022 0.639809 0.54324 0.792887 0.547553 0.839619 0.200446 0.476359 0.427033 0.210896 0.37966 0.194706 0.439351 0.594902 ~~~~ Julia s0 = RandomRuta(D) coste_s0 = coste(s0) ~~~~ Julia output Ruta([39, 49, 8, 12, 40, 10, 22, 46, 43, 25 … 35, 20, 27, 6, 13, 17, 34, 1, 26, 24], [0.0 0.3709524487108521 … 0.7061405535731123 0.3412001689877821; 0.3709524487108521 0.0 … 0.38490283213267174 0.6892542891509004; … ; 0.7061405535731123 0.38490283213267174 … 0.0 0.9600886854862052; 0.3412001689877821 0.6892542891509004 … 0.9600886854862052 0.0]) 30.047661586127536 ~~~~ Julia # Ejecuta Templado Simulado sol_ruta, sol_coste = SA(coste, s0, 1000, cambia_tour, 0.0001, .1, true) ~~~~ Julia output (Ruta([41, 49, 36, 28, 38, 27, 22, 29, 47, 24 … 7, 35, 14, 17, 21, 33, 32, 6, 40, 4], [0.0 0.3709524487108521 … 0.7061405535731123 0.3412001689877821; 0.3709524487108521 0.0 … 0.38490283213267174 0.6892542891509004; … ; 0.7061405535731123 0.38490283213267174 … 0.0 0.9600886854862052; 0.3412001689877821 0.6892542891509004 … 0.9600886854862052 0.0]), 6.095636154696887) ~~~~ Julia # Mostramos el resultado gráficamente # Reordenamos las ciudades según el orden encontrado en la solución ciu_orden = map(o -> ciudades[:,o], sol_ruta.orden) # Copiamos la última ciudad y la ponemos en primer lugar ciu_orden = pushfirst!(ciu_orden,ciu_orden[num_ciudades]) # Añadimos las ciudades originales como puntos scatter(ciudades[1,:], ciudades[2,:], color=:red, markersize=5) # Las numeramos según el orden inicial annotate!(ciudades[1,:], ciudades[2,:], [text(i,:bottom) for i in 1:num_ciudades]) # Dibujamos la ruta solución plot!(map(c -> c[1], ciu_orden),map(c -> c[2], ciu_orden),color=:blue) # Devolvemos la figura ~~~~ 

## Mejoras al Templado Simulado Hay algunas mejoras inmediatas que pueden hacerse sobre este algoritmo si tenemos en cuenta casos particulares de su aplicación. Debe tenerse en cuenta que, por ejemplo, el algoritmo así planteado no restringe su utilidad al caso de un número finito ni discreto de acciones, y se puede usar tanto para búsquedas continuas como discretas. ### Caso continuo Además de utilizarse en la búsqueda en espacios discretos (**optimización combinatoria**), el templado simulado puede extenderse a la minimización de funciones continuas, donde si hay mucha variabilidad es muy fácil quedarse atascado en uno de los muchos mínimos locales. El problema empeora a medida que la dimensionalidad del problema (el número de variables que intervienen) crece. El Templado Simulado, comparado con otros métodos habituales en búsqueda optimizada continua (como el [método simplex de Nelder-Mead](https://es.wikipedia.org/wiki/M%C3%A9todo_Nelder-Mead), o la [búsqueda aleatoria adaptativa GRASP](https://ccc.inaoep.mx/~emorales/Cursos/Busqueda/node66.html)) sobre una serie de funciones continuas de prueba, da muy buenos resultados, consiguiendo identificar el mínimo global, y de forma más eficiente si comparamos el número de evaluaciones de la función que son necesarias, pero sigue siendo un procedimiento costoso desde el punto de vista computacional, ya que el número de evaluaciones crece linealmente con el número de dimensiones. A pesar de ello, este coste es menos prohibitivo que el de las otras técnicas, sobre todo porque el coste depende únicamente de la temperatura inicial, lo que significa que el tiempo de cálculo se conoce de antemano, una ventaja del algoritmo. !!!code:Julia

Para hacer una optimización de una función continua basta considerar pequeñas variaciones del dato de entrada. Lo más probbale es que no devuelva el resultado exacto, pero se pueden conseguir buenas aproximaciones. Por ejemplo, imaginemos que queremos calcular los mínimos de una función real, $f:\mathbb{R} \to \mathbb{R}$, dada por: $$f(x) = 1+ cos(s)^2 - sin(x)$$ y hagamos dos intentos, uno partiendo de $x= 0$ y otro de $x=2.5$ ~~~~Julia f(x) = 1 + cos(x)^2 - sin(x) vecino(x) = x + rand([-0.01, 0.01]) x1, y1 = SA(f, 0, 100, vecino, 0.001, 1.0, true) ~~~~ Julia output (1.5700000000000012, 9.512043955828275e-7) ~~~~Julia x2, y2 = SA(f, -2.5, 100, vecino, 0.001, 1.0, false) ~~~~ Julia output (-4.709999999999953, 8.56082870337449e-6) ~~~~ Podemos dar una representación gráfica de $f$ en el intervalo $[-2\pi,2\pi]$ con los puntos calculados superpuestos, para comprobar el resultado: ~~~~Julia In = -2π:0.1:2π plot(In, map(f, In), label="f(x)") plot!([x1,x2], [y1,y2], seriestype=:scatter, color = [:red, :green], markersize = 5, label = "SA") ~~~~ 

### Paralelización Uno de los principales problemas del Templado Simulado, tanto para el caso discreto como continuo, es el elevado tiempo de cálculo para encontrar una solución cercana a un óptimo. Este problema surge al proponer y evaluar los valores de las funciones para luego rechazarlos y puede mitigarse en cierta medida paralelizando el proceso. Muchos intentos de paralelización del Templado Simulado para problemas específicos han tenido éxito, especialmente para problemas bien estudiados (como el problema del viajante de negocios). Esta paralelización suele hacerse dividiendo el problema en subproblemas más pequeños, cada uno de los cuales se resuelve mediante Templado Simulado estándar, como el presentado aquí, y luego recombinando estas soluciones. Dividir el problema de esta manera debe hacerse con mucho cuidado, ya que la alta interdependencia de los subproblemas puede ser problemática. Si los subproblemas interdependientes no se comunican, la recombinación conducirá a una solución que está potencialmente lejos del óptimo global, análoga a las estructuras localmente óptimas del enfriamiento del sistema en su conjunto. Por otra parte, si los subproblemas se comunican, debido a su gran interdependencia, se requiere la comunicación entre los procesadores, lo que sólo supone pequeñas ganancias de eficiencia en comparación con la ejecución de un templado secuencial. !!!note: Clases de Problemas y Paralelización Los métodos de paralelización diferencian dos clases de problemas para medir si compensa paralelizarlos o no cuando se aplica el Templado Simulado: * Clase $K_1$: el tiempo para generar configuraciones vecinas es aproximadamente igual al tiempo que se tarda en calcular el incremento de energía entre la configuración actual y la propuesta. Por ejemplo: el problema de viajante de negocios. * Clase $K_2$: el tiempo para generar una configuración vecina es mucho menor que el tiempo para calcular el incremento de energía. Por ejemplo: el problema clustering (que veremos en la parte de aprendizaje automático). # Optimización por Enjambres de Partículas La **Optimización por Enjambres de Partículas** (conocida como **PSO**, por sus siglas en inglés, **Particle Swarm Optimization**) es otro algoritmo de optimización que se diferencia del Templado Simulado en dos puntos destacables: * el espacio de estados se puede representar como un espacio vectorial $n$-dimensional, $\mathbb{R}^n$, y * usa un conjunto de buscadores (que aquí se llaman *partículas*) para explorar el espacio de estados siguiendo trayectorias muy simples y que ponen en común sus hallazgos (podríamos entender que en Templado Simulado solo hay un buscador que, además, se mueve de forma estocástica por el espacio). !!!side:4 Kennedy, J., Eberhart, R. _Particle swarm optimization_, Neural Networks, 1995. Proceedings IEEE International Conference. ![](img/abeja.jpg align="right" width=20%)Este método fue descrito alrededor de 1995 por Kennedy y Eberhart $^4$, y de hecho se inspiraron en el comportamiento de los enjambres de insectos en la naturaleza para su definición formal. En concreto, podemos imaginar un enjambre de abejas que, a la hora de buscar polen, recorren la región del espacio en la que existe más densidad de flores, porque la probabilidad de que haya polen es mayor, y a la vuelta a la colmena comunican al resto de individuos los resultados de su recorrido. La misma idea fue trasladada al campo de la computación en forma de algoritmo y se emplea en la actualidad en la optimización de distintos tipos de sistemas. Como en el resto de problemas de optimización, se supone que tenemos una función desconocida, $f(\mathbf{x})$ (fitness u objetivo), que podemos evaluar en los puntos que queramos pero a modo de caja negra (sin disponer de una expresión analítica), y el objetivo es encontrar valores de $\mathbf{x}$ para los que la función $f(\mathbf{x})$ tenga un óptimo (máximo, mínimo, o bien verifica alguna relación extremal respecto a alguna otra función). La idea que vamos a seguir en PSO comienza situando partículas al azar en el espacio de búsqueda, pero dándoles la posibilidad de que se muevan a través de él de acuerdo a unas reglas deterministas y simples que tienen en cuenta el conocimiento personal de cada partícula y el conocimiento global del enjambre acerca del espacio de estado. Veremos que proporcionándoles una capacidad simple de movimiento por este espacio y de comunicación entre ellas pueden llegar a descubrir valores particularmente buenos para $f(\mathbf{x})$ gastando pocos recursos computacionales (cálculos, memoria y tiempo). Para fijar los elementos que usaremos, supondremos que tenemos $N$ partículas, denotadas $\{1,\dots,N\}$, y cada partícula tiene una **posición**, $\mathbf{x}_i$, en el espacio de búsqueda y una **velocidad**, $\mathbf{v}_i$, que determina su movimiento a través del espacio. Además, simulando partículas de un mundo real físico (y de aquí viene el hecho de llamarlas **partículas**), tienen una cantidad de **inercia** que las mantiene en la misma dirección en la que se movían, y una **aceleración** (cambio de velocidad) que depende principalmente de dos características: 1. Cada partícula es atraída hacia la mejor localización que ella, personalmente, ha encontrado en su historia (**mejor personal**), $\mathbf{x}_i^{pb}$. !!!side:5 Como veremos más adelante, y en la implementación en Julia del algoritmo, no es necesario que cada partícula tenga comunicación con el enjambre completo, sino con un sbconjunto de ellas que actúan a modo de _vecindario_ social. 2. Cada partícula es atraída hacia la mejor localización que ha sido encontrada por el conjunto de partículas $^5$ en el espacio de búsqueda (**mejor global**), $\mathbf{x}^{gb}$.  La fuerza con que las partículas son empujadas en cada una de estas direcciones depende de dos parámetros que pueden ajustarse (**atracción-al-mejor-personal** y **atracción-al-mejor-global**) de foma que, a medida que las partículas se alejan de estas localizaciones mejores, la fuerza de atracción es mayor. También se suele incluir un factor aleatorio que influye en cómo las partículas son empujadas hacia estas localizaciones. Formalmente, lo podemos escribir como: $$\begin{equation} \mathbf{v}_i(t+1)=\mathbf{v}_i(t)+c_1\cdot r_1\cdot (\mathbf{x}_i^{pb}−\mathbf{x}_i(t))+c_2\cdot r_2\cdot (\mathbf{x}^{gb}−\mathbf{x}_i(t)) \end{equation} $$ donde: - las posiciones y velocidades se indican en función del tiempo (discreto, por pasos), - $c_1$ y $c_2$ son, respectivamente, las constantes de atracción al mejor personal y mejor global, y - $r_1$ y $r_2$ son dos números aleatorios en $[0,1]$ (que se lanzan en cada iteración y para cada partícula). Una vez actualizadas las velocidades de todas las partículas, sus posiciones se actualizan siguiendo una ley simple de movimiento: $$ \begin{equation} \mathbf{x}_i(t+1)=\mathbf{x}_i(t)+\mathbf{v}_i(t) \end{equation} $$  El algoritmo PSO se esboza a grandes rasgos en el siguente código: ```none Asignar posiciones y velocidades aleatorias iniciales a las partículas Repetir Para cada partícula del enjambre: Actualizar v haciendo uso de (1) Actualizar x haciendo uso de (2) Calcular f(x) Actualizar mejor personal Actualizar mejor global Devolver x^gb ```  Lo interesante de este modelo es que cada partícula solo hace operaciones vectoriales básicas (para actualizar las velocidades y posiciones basta realizar algunas sumas de vectores) y, como la experiencia muestra que hacen falta pocas partículas y el número de iteraciones es bajo, la complejidad del proceso se centra sobre todo en la evaluación de $f$ de cada partícula en cada iteración, algo que no puede evitarse y que depende exclusivamente del problema específico que se quiera optimizar. En consecuencia, este método de optimización es especialmente adecuado cuando el coste de la función de fitness es muy elevado, pues requiere pocas evaluaciones para encontrar valores cercanos al óptimo. Como inconveniente, no podemos estar seguros de alcanzar el óptimo global, como pasa en todo este tipo de metaheurísticas (problema que comparten todos los métodos de optimización que trabajan con cajas negras). !!!side:8 El valor de la función no hay necesidad de modificarlo. !!!side:9 Si $\pi$ es una transformación lineal (habitualmente será un proceso de escala y traslación en cada coordenada), entonces $\pi^{-1}$ también es una transformación lineal que se puede calcular fácilmente a partir de $\pi$. !!!note:Estandarizando la implementación Respecto a una posible implementación, es habitual estandarizar los valores posibles que puede tomar tanto la posición de las partículas como sus velocidades. Concretamente, es habitual que $\mathbf{x}_i\in [0,1]^n$, mientras que $\mathbf{v}_i\in [-1,1]^n$. En caso de que la función no esté definida en este espacio (que será lo más habitual) se procede a hacer una transformación entre el espacio de definición del problema original y el que admite el algoritmo, y posteriormente se aplica la transformación inversa a la posición de la mejor partícula encontrada $^8$. Si el problema está definido en $D\subseteq \mathbb{R}^n$, se considera una transformación: $$\pi:[0,1]^n \to D$$ y se aplica PSO a la función $f'=f \circ \pi$. Tras encontrar $\mathbf{x}^{gb}$ para $f'$, se devuelve $\pi^{-1}(\mathbf{x}^{gb})$ como mejor aproximación $^9$ para $f$. Además de facilitar la implementación, tener todas las coordenadas moviéndose en intervalos de igual tamaño hace que la búsqueda sea más homogénea, ya que la diversidad de escalas puede provocar que haya mayor tendencia a moverse en direcciones en las que el tamaño del parámetro correspondiente puede ser mayor. !!!code:Julia

Vamos a definir todas las estructuras y funciones auxiliares que son necesarias para poder aplicar el modelado y solución por medio de PSO. En general, recordemos que estamos intentando minimizar una función, $f$, que podemos evaluar puntualmente (pero no es necesario conocer su expresión). Por tanto, el primer paso es definir la estructura que contendrá la información cada partícula: - `nDim::Int`: dimensión del espacio de parámetros a explorar - `pos::Array{Float64, 1}`: posición actual de la partícula - `vel::Array{Float64, 1}`: velocidad actual de la partícula - `pBest::Array{Float64, 1}`: posición en la que la partícula tiene el valor que mejor se ajusta a través de la historia de la partícula - `gBest::Array{Float64, 1}`: posición en la que el grupo local tiene el valor que mejor se ajusta a través de la historia del grupo local - `fitValue::Float64`: valor de f en la posición actual - `fitpBest::Float64`: valor de f en `pBest` - `fitgBest::Float64`: valor de f en `gBest` Observa que en la implementación se proporciona también un constructor personalizado para la estructura, que permite declarar una partícula únicamente usando su dimensión, ya que el resto de valores son aleatorios: ~~~~ Julia linenumbers mutable struct Particle nDim::Int pos::Array{Float64, 1} vel::Array{Float64, 1} pBest::Array{Float64, 1} gBest::Array{Float64, 1} fitValue::Float64 fitpBest::Float64 fitgBest::Float64 # Constructor function Particle(nDim::Int) pos = rand(nDim) # Posición aleatoria (∈ [0,1]^n) vel = rand(nDim) - pos # Velocidad aleatoria (∈ [-1,1]^n) pBest = gBest = pos # Memoria inicial fitValue = fitpBest = fitgBest = Inf # fitValue = fitpBest = fitgBest = ∞ new(nDim, pos, vel, pBest, gBest, fitValue, fitpBest, fitgBest) end end # p = Particle(4) ~~~~ Las funciones de actualización de la partícula son implementaciones directas de algoritmo mostrado más arriba: ~~~~ Julia linenumbers # Actualiza `fitValue` de `p` utilizando la función `fitFunc`. function initFitValue!(fitFunc::Function, p::Particle) p.fitValue = fitFunc(p.pos) nothing end # Actualiza posición y `fitValue` de la partícula `p` en la nueva posición. function updatePosAndFitValue!(fitFunc::Function, p::Particle) p.pos += p.vel # si la posición está fuera del espacio de parámetros, [0,1]^n, # establecemos fitValue = Inf for x in p.pos if (x < 0 || x > 1) p.fitValue = Inf return end end # actualiza valor p.fitValue = fitFunc(p.pos) nothing end # Actualiza `pBest` y `fitpBest` para la partícula `p`. function updatepBestAndFitpBest!(p::Particle) if p.fitValue < p.fitpBest p.fitpBest = p.fitValue p.pBest = p.pos end nothing end # Actualiza `vel` de la partícula `p`. function updateVel!(p::Particle, w::Float64, c1::Float64, c2::Float64) p.vel = w * p.vel + c1 * rand() * (p.pBest - p.pos) + c2 * rand() * (p.gBest - p.pos) nothing end ~~~~ Para que la implementación sea un poco más general que la mostrada, vamos a permitir que los vecinos de cada partícula se definan en forma un poco más particular. En concreto, supondremos que los índices de las partículas forman una estructura cíclica: `1,2,...,n,1,2...`, y los vecinos de cada partícula se definen como un radio alrededor de su índice ~~~~ Julia linenumbers # Devuelve los índices de los `nNeighbor` vecinos de la partícula `i` en un enjambre con `nParticle` partículas tomando `nNeighbor` vecinos en cada dirección. function neighIndices(i::Int, nNeighbor::Int, nParticle::Int) # el número de vecinos debe ser superior a 3 nNeighbor = max(3, nNeighbor) # número de vecinos a la izquierda de la partícula i-ésima nLeft = (nNeighbor - 1) ÷ 2 # el índice de la primera partícula del grupo local startIndex = (i - nLeft) # el índice de la última partícula del grupo local endIndex = startIndex + nNeighbor -1 # índices para el grupo local indices = collect(startIndex:endIndex) # ajusta los índices para que estén en range(1:nParticle) for i in 1:nNeighbor if indices[i] < 1 indices[i] += nParticle elseif indices[i] > nParticle indices[i] -= nParticle end end indices end # neighIndices(1, 5, 40) ~~~~ Definamos ahora la estructura que declara un enjambre, y que contiene los siguientes campos (de nuevo, damos un constructor personalizado): - `fitFunc::Function`: función que debe evaluarse - `nDim::Int`: dimensión del espacio de parámetros a explorar - `nParticle::Int`: número de partículas de un enjambre - `nNeighbor::Int`: número de vecinos (partículas) en un grupo local - `nIter::Int`: número de iteraciones de actualización del ejambre - `c1::Float64`: constante cognitiva - `c2::Float64`: constante social - `wMax::Float64`: valor máximo del peso de la inercia - `wMin::Float64`: valor mínimo del peso de la inercia - `w::Float64`: valor actual del peso de inercia - `gBest::Array{Float64, 1}`: posición en la que el enjambre tiene el valor mejor a lo largo de la historia - `fitgBest::Float64`: valor en `gBest` - `particles::Array{Particle, 1}`: partículas del enjambre ~~~~ Julia linenumbers mutable struct Swarm fitFunc::Function nDim::Int nParticle::Int nNeighbor::Int nIter::Int c1::Float64 c2::Float64 wMax::Float64 wMin::Float64 w::Float64 gBest::Array{Float64, 1} fitgBest::Float64 particles::Array{Particle, 1} # Constructor function Swarm(fitFunc::Function, nDim::Int; nParticle::Int=40, nNeighbor::Int=3, nIter::Int=2000, c1::Float64=2.0, c2::Float64=2.0, wMax::Float64=0.9, wMin::Float64=0.4) # El tamaño del vecindario no puede superar el del enjambre if nNeighbor > nParticle nNeighbor = nParticle end w = wMax gBest = rand(nDim) fitgBest = Inf # Crear el enjambre con nPartículas particles = [Particle(nDim) for i in 1:nParticle] new(fitFunc, nDim, nParticle, nNeighbor, nIter, c1, c2, wMax, wMin, w, gBest, fitgBest, particles) end end # f(x)=x^2+1 # s = Swarm(f, 2) ~~~~ Y definimos las funciones que permiten actualizar un paso de ejecución del algoritmo el enjambre, así inicializar el enjambre: ~~~~ Julia linenumbers # Actualiza gBest y fitgBest para cada partícula del enjambre `s`. function updategBestAndFitgBestParticles!(s::Swarm) for i in 1:s.nParticle neighborIds = neighIndices(i, s.nNeighbor, s.nParticle) neighborFits = [s.particles[Id].fitValue for Id in neighborIds] fitgBest, index = findmin(neighborFits) if fitgBest < s.particles[i].fitgBest gBest = s.particles[neighborIds[index]].pos s.particles[i].gBest = gBest s.particles[i].fitgBest = fitgBest end end nothing end # Actualiza gBest y fitgBest para el enjambre `s`. function updategBestAndFitgBest!(s::Swarm) gFits = [particle.fitValue for particle in s.particles] fitgBest, index = findmin(gFits) if fitgBest < s.fitgBest s.gBest = s.particles[index].pos s.fitgBest = fitgBest end nothing end # Inicialización (0ª iteración) del enjambre `s`. function initSwarm(s::Swarm) # Reiniciar el fitValue para cada partícula for particle in s.particles initFitValue!(s.fitFunc, particle) updatepBestAndFitpBest!(particle) end # actualizar gBest y fitgBest para las partículas del enjambre updategBestAndFitgBestParticles!(s) # actualizar gBest y fitgBest para el enjambre updategBestAndFitgBest!(s) nothing end # Actualiza el peso de inercia después de cada iteración. function updateInertia!(s::Swarm) dw = (s.wMax - s.wMin)/s.nIter s.w -= dw nothing end # Actualiza las partículas del enjambre `s`. function updateParticles!(s::Swarm) for particle in s.particles updateVel!(particle, s.w, s.c1, s.c2) updatePosAndFitValue!(s.fitFunc, particle) end nothing end ~~~~ A partir de las funciones anteriores, un paso de iteración del enjambre `s` sería: ~~~~ Julia linenumbers function updateSwarm!(s::Swarm) # actualiza velocidad, posición y fitness de las partículas del enjambre updateParticles!(s::Swarm) # actualizar los valores gBest y fitgBest de cada partícula del enjambre updategBestAndFitgBestParticles!(s::Swarm) # actualizar el gBest y fitgBest para el enjambre updategBestAndFitgBest!(s::Swarm) # actualizar el peso de inercia w para cada partícula del enjambre updateInertia!(s::Swarm) nothing end ~~~~ La función que permite crear un enjambre con las características pedidas y ejecutar PSO con ese enjambre `nIter` pasos es ya sencilla usando la librería de funciones que hemos definido. Esta función dvuelve la mejor posición encontrada, y el valor de la función de optimización en esa posición. ~~~~ Julia linenumbers function swarmFit(fitFunc::Function, nDim::Int; nParticle::Int = 10, nIter::Int = 100, nNeighbor::Int = 3) # Crea el ejambre y lo inicializa s = Swarm(fitFunc, nDim, nParticle=nParticle, nIter=nIter, nNeighbor=nNeighbor) initSwarm(s) # Ejecuta la nIter actualizaciones del enjambre for i in 1:s.nIter updateSwarm!(s) end # Devuelve la posición de la mejor solución encontrada return s.gBest, s.fitgBest end ~~~~ Vamos a mostrar un ejemplo de uso muy sencilla. Supongamos que queremos calcular los mínimos de la función 2D: $$ f(x,y) = x^2 + y^2$$ Definamos para ello la función en Julia, y un enjambre de 4 partículas de dimensión 2: ~~~~ Julia # Observa que definimos f como una función sobre vectores 2D f(x) = x[1]^2 + x[2]^2 s = Swarm(f, 2, nParticle=4) ~~~~ Julia output Swarm(f, 2, 4, 3, 2000, 2.0, 2.0, 0.9, 0.4, 0.9, [0.46403539901495605, 0.24500219608925766], Inf, Particle[Particle(2, [0.18835787094253886, 0.7222892509090975], [0.70896380285857, 0.2682851420309852], [0.18835787094253886, 0.7222892509090975], [0.18835787094253886, 0.7222892509090975], Inf, Inf, Inf), Particle(2, [0.04054666451861566, 0.7136157307850581], [0.3003446933100331, 0.08136599096328245], [0.04054666451861566, 0.7136157307850581], [0.04054666451861566, 0.7136157307850581], Inf, Inf, Inf), Particle(2, [0.23138771768402067, 0.370212650032125], [0.026761377468961123, 0.4206467964867012], [0.23138771768402067, 0.370212650032125], [0.23138771768402067, 0.370212650032125], Inf, Inf, Inf), Particle(2, [0.8499730245086784, 0.9688459080739353], [-0.7964170642294439, 0.011056752737428233], [0.8499730245086784, 0.9688459080739353], [0.8499730245086784, 0.9688459080739353], Inf, Inf, Inf)]) ~~~~ Julia # Inicializamos el enjambre initSwarm(s) # y mostramos la primera partícula del enjambre (que debe tener datos aleatorios) println(s.particles[1]) ~~~~ Julia output Particle(2, [0.18835787094253886, 0.7222892509090975], [0.70896380285857, 0.2682851420309852], [0.18835787094253886, 0.7222892509090975], [0.04054666451861566, 0.7136157307850581], 0.5571804495248313, 0.5571804495248313, 0.5108914432274778) ~~~~ Julia # Actualicemos una vez el enjambre para comprobar que, en efecto, la partícula se ha movido: updateSwarm!(s) println(s.particles[1]) ~~~~ Julia output Particle(2, [0.5499373666310641, 0.9475216447570105], [0.36157949568852527, 0.22523239384791302], [0.18835787094253886, 0.7222892509090975], [0.4463797184860497, 0.5429827221964985], 1.2002283745001399, 0.5571804495248313, 0.4940850896796048) ~~~~ Julia # Ejecutemos ya el algoritmo PSO para un enjambre adecuado: con 100 partículas, cada una mirando a 10 vecinos a cada lado, y dando 50 pasos: # swarmFit(f, 2, nParticle = 100, nIter = 50, nNeighbor = 10) ~~~~ Julia output ([0.015536459142767467, 0.011828390915045062], 0.00038129239433400337) ~~~~ Tras $50$ iteraciones, el enjambre encuentra un mínimo en $(0.01, 0.01)$ y con un valor de $f$ de $0.0003$. El mínimo real de $f$ es $0$ y se encuentra en el $(0,0)$.

!!!warn: ¡¡¡Atención!!! Recuerda que la librería está preparada para trabajar con partículas definidas en $[0,1]^D$, por lo que si tu problema necesita partículas en otros dominios, debes modificarlo previamente (haciendo un escalado lineal) para que se adapte a la implementación. Y no olvides realizar el escalado inverso sobre el resultado para poder saber qué posición ocupa el mínimo en el dominio de tu problema. ## Variantes Una primera posibilidad, que se implementa en algunos de los algoritmos que se engloban dentro de los PSO, y que es fácil de implementar en el modelo anterior, consiste en añadir una fuerza repulsiva entre las partículas para intentar prevenir que todas ellas converjan prematuramente a una pequeña zona del espacio de búsqueda (lo que supone acabar prematuramente en un óptimo local, y con pocas opciones de escapar de él para encontrar mejores optimizaciones). Se puede enriquecer el modelo añadiendo comportamientos diversos a los individuos-partículas, ya sea inspirándose en comportamientos observados en la naturaleza o puramente abstractos, y hay diversos ejemplos de ello en la literatura relacionada con este tipo de algoritmos de optimización. Como ya comentamos en una nota anterior, y se implementa en cierta forma en la versión dada en Julia, otra opción es no considerar el óptimo global que se va encontrando, sino dividir las partículas en familias/vecindarios y considerar óptimos globales para estos grupos. De esta forma, cada familia puede trabajar en paralelo y cada partícula solo ha de tener conocimiento acerca de lo que *sucede* en su familia, no en el conjunto global de individuos. !!!side:10 Una variante habitual de este último tipo es considerar que el vecindario de una partícula viene dada por las partículas que están dentro de un determinado entorno (por ejemplo, a menos de una determinada distancia). Además, existen variantes en las que estos vecindarios son dinámicos, pueden ir cambiando a medida que el algoritmo se ejecuta $^{10}$. Así como permitir que cada partícula pertenezca a más de un vecindario, lo que permite que haya comunicaciones a larga distancia entre partículas aunque no pertenezcan al mismo vecindario (como es el caso de la implementación que hemos dado). Por último, ¿qué pasa si la función que intenta ser optimizada cambia en el tiempo?, es decir, en el caso en que sea una búsqueda dinámica, en la que la variación de la función durante el tiempo de ejecución haga que los extremos se vayan desplazando por el espacio de búsqueda. Normalmente, los enjambres funcionan bien incluso bajo el funcionamiento de lo que se denominan **paisajes dinámicos**.