Optimización con enjambre de partículas (Particle Swarm Optimization)
Más sobre ciencia de datos: cienciadedatos.net
Introducción¶
La optimización por enjambre de partículas (Particle Swarm Optimization, PSO) es un método de optimización heurística orientado a encontrar mínimos o máximos globales. Su funcionamiento está inspirado en el comportamiento que tienen las bandadas de pájaros o bancos de peces en los que, el movimiento de cada individuo (dirección, velocidad, aceleración...), es el resultado de combinar las decisiones individuales de cada uno con el comportamiento del resto.
El método de enjambre de partículas es solo una de las muchas estrategias de optimización heurística que existen, una alternativa común son los algoritmos genéticos.
La optimización heurística no tiene por qué ser la forma de optimización más adecuada en todos los escenarios. Si el problema en cuestión puede optimizarse de forma analítica, suele ser más adecuado resolverlo de esta forma.
La implementación de algoritmo que se muestra en este documento pretende ser lo más explicativa posible aunque para ello no sea la más eficiente.
El código de las funciones desarrolladas a lo largo del documento puede descargarse en el siguiente link.
Algoritmo¶
Aunque existen variaciones, algunas de las cuales se describen a lo largo de este documento, en términos generales, la estructura de un algoritmo PSO para optimizar (maximizar o minimizar) una función con una o múltiples variables sigue los siguientes pasos:
Crear un enjambre inicial de $n$ partículas aleatorias. Cada partícula consta de 4 elementos: una posición que representa una determinada combinación de valores de las variables, el valor de la función objetivo en la posición donde se encuentra la partícula, una velocidad que indica cómo y hacia donde se desplaza la partícula, y un registro de la mejor posición en la que ha estado la partícula hasta el momento.
Evaluar cada partícula con la función objetivo.
Actualizar la posición y velocidad de cada partícula. Esta es la parte que proporciona al algoritmo la capacidad de optimización. En el apartado Mover partícula se describe con detalle el proceso.
Si no se cumple un criterio de parada, volver al paso 2.
En los siguientes apartados se implementan cada una de las etapas del proceso para, finalmente, combinarlas todas en una única función.
Crear partícula¶
Cada partícula está definida por una posición, velocidad y valor que varían a medida que la partícula se mueve. Además, también almacena la mejor posición en la que ha estado hasta el momento. Cuando se crea aun nueva partícula, únicamente se dispone de información sobre su posición y velocidad (normalmente iniciada como cero), el resto de valores no se conocen hasta que la partícula es evaluada.
Evaluar partícula¶
Evaluar una partícula consiste en calcular el valor de la función objetivo en la posición que ocupa la partícula es ese momento. Cada partícula almacena también la posición con mejor valor en la que ha estado hasta el momento. Para poder identificar si una nueva posición es mejor que las anteriores, es necesario conocer si se trata de un problema de minimización o maximización.
Mover partícula¶
Mover una partícula implica actualizar su velocidad y posición. Este paso es el más importante ya que otorga al algoritmo la capacidad de optimizar.
La velocidad de cada partícula del enjambre se actualiza empleando la siguiente ecuación:
$$v_i(t+1) = wv_i(t) + c_1r_1[\hat{x}_i(t) - x_i(t)] + c_2r_2[g(t) - x_i(t)]$$donde:
- $v_i(t+1)$: velocidad de la partícula $i$ en el momento $t + 1$, es decir, la nueva velocidad.
- $v_i(t)$: velocidad de la partícula $i$ en el momento $t$, es decir, la velocidad actual.
- $w$: coeficiente de inercia, reduce o aumenta a la velocidad de la partícula.
- $c_1$: coeficiente cognitivo.
- $r_1$: vector de valores aleatorios entre 0 y 1 de longitud igual a la del vector velocidad.
- $\hat{x}_i(t)$: mejor posición en la que ha estado la partícula $i$ hasta el momento.
- $x_i(t)$: posición de la partícula $i$ en el momento $t$.
- $c_2$: coeficiente social.
- $r_2$: vector de valores aleatorios entre 0 y 1 de longitud igual a la del vector velocidad.
- $g(t)$: posición de todo el enjambre en el momento $t$, el mejor valor global.
Para comprender como se relaciona esta ecuación con el movimiento de la partícula, resulta útil diferenciar tres partes:
$wv_i(t)$ es la componente de inercia, responsable de mantener a la partícula moviéndose en la dirección en la que lo ha estado haciendo hasta el momento. El valor recomendado del coeficiente de inercia $w$ suele ser entre 0.8 y 1.2. Si $w<1$, la partícula se va desacelerando a medida que avanzan las iteraciones, esto se traduce en menor exploración pero una convergencia hacia el óptimo más rápida. Si $w>1$, la partícula se va acelerando, lo que permite explorar más zonas del espacio de la función pero dificulta la convergencia.
$c_1r_1[\hat{x}_i(t) - x_i(t)]$ es la componente cognitiva, responsable de que la partícula tienda a moverse hacia la posición donde ha obtenido mejores resultados hasta el momento. El coeficiente cognitivo $c_1$ suele estar acotado en el rango [0, 2], siendo 2 el valor recomendado. $r_1$ es un vector de valores aleatorios entre 0 y 1 (un valor por cada dimensión) que aporta cierto comportamiento estocástico al movimiento de las partículas, mejorando así la capacidad de escapar de mínimos locales.
$c_2r_2[g(t) - x_i(t)]$ es la componente social, responsable de que la partícula tienda a moverse hacia la mejor posición encontrada por el enjambre hasta el momento. Puede interpretarse como el "conocimiento colectivo". El valor del coeficiente social $c_2$ suele estar acotado en el rango [0, 2], siendo 2 el valor recomendado. $r_2$ es un vector de valores aleatorios entre 0 y 1 (un valor por cada dimensión) que aporta cierto comportamiento estocástico al movimiento de las partículas, mejorando así la capacidad de escapar de mínimos locales.
La magnitud relativa entre la componente cognitiva y la componente social permite regular el comportamiento exploratorio del algoritmo. Cuanto mayor es el valor de $c_1$ respecto a $c_2$, mayor independencia de movimiento tiene cada partícula, lo que permite mayor exploración pero mayor lentitud en la convergencia. Por el contrario, cuanto mayor es el valor de $c_2$ respecto a $c_1$, más obligadas están las partículas a moverse hacia la mejor región encontrada hasta el momento, lo que reduce la exploración pero acelera la convergencia.
En algunas versiones del algoritmo, $r_1$ y $r_2$ son escalares en lugar de vectores. Multiplicar cada componente de la velocidad por un valor aleatorio distinto añade mayores fluctuaciones al movimiento de las partículas, lo que, aun a riesgo de retrasar la convergencia, suele generar mejores resultados.
Una vez calculada la nueva velocidad, se puede actualizar la posición de la partícula con la ecuación:
$$x_i(t+1) = x_i(t) + v_i(t+1)$$Uno de los principales problemas del algoritmo PSO es que las partículas suelen adquirir velocidades excesivamente altas, lo que les lleva a salirse de los límites del espacio de búsqueda o a que sean incapaces de converger en la región óptima. Es en este paso del algoritmo donde más investigaciones y adaptaciones se han hecho. Algunas de las soluciones son:
Limitar la velocidad máxima que puede alcanzar una partícula. Siendo [$x_{min}$, $x_{max}$] los límites inferior y superior del espacio de búsqueda de cada variable, la velocidad máxima que puede alcanzar la partícula en esa dirección es $v_{max} = k(x_{max} - x_{min})/2$, donde $k$ suele ser un valor entre 0.1 y 1.
Si el valor de alguna de las variables excede los límites impuestos, se sobrescribe con el valor del límite correspondiente y se reinicia su velocidad a cero.
Reducción lineal del coeficiente de inercia $w$. Esta estrategia consiste en ir reduciendo el coeficiente de inercia a medida que avanzan las iteraciones. En las primeras iteraciones, las partículas tiene mucha capacidad de exploración y, a medida que avanza el proceso, va reduciéndose su velocidad favoreciendo la convergencia. Puede conseguirse este efecto con la ecuación:
donde:
- $w_{t}$: coeficiente de inercia en la iteración $t$.
- $w_{max}$: coeficiente de inercia máximo. Valor con el que se inicia el algoritmo. Valor recomendado de 0.9.
- $w_{min}$: coeficiente de inercia mínimo. Valor que se alcanza en la última iteración. Valor recomendado de 0.4.
- $t_{max}$: número máximo de iteraciones.
La siguiente función actualiza la posición de una partícula teniendo en cuenta su posición y velocidad actual, la mejor posición global encontrada por el enjambre, los coeficientes de inercia, cognitivo, social, y los límites de búsqueda.
Implementación python¶
################################################################################
# LIBRERÍAS NECESARIAS
################################################################################
import numpy as np
import random
import warnings
import random
import copy
import pandas as pd
import time
import matplotlib
import matplotlib.pyplot as plt
from datetime import datetime
# Configuración warnings
################################################################################
import warnings
warnings.filterwarnings('ignore')
Clase partícula¶
################################################################################
# CLASE PARTÍCULA #
################################################################################
class Particula:
"""Esta clase representa nueva partícula con una posición inicial definida por
una combinación de valores numéricos aleatorios y velocidad de 0. El rango
de posibles valores para cada variable (posición) puede estar acotado. Al
crear una nueva partícula, solo se dispone de información sobre su posición
inicial y velocidad, el resto de atributos están vacíos.
Parameters
----------
n_variables : `int`
número de variables que definen la posición de la partícula.
limites_inf : `list` or `numpy.ndarray`, optional
límite inferior de cada variable. Si solo se quiere predefinir límites
de alguna variable, emplear ``None``. Los ``None`` serán remplazados
por el valor (-10**3). (default is ``None``)
limites_sup : `list` or `numpy.ndarray`, optional
límite superior de cada variable. Si solo se quiere predefinir límites
de alguna variable, emplear ``None``. Los ``None`` serán remplazados
por el valor (+10**3). (default is ``None``)
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
Attributes
----------
n_variables : `int`
número de variables que definen la posición de la partícula.
limites_inf : `list` or `numpy.ndarray`
límite inferior de cada variable. Si solo se quiere predefinir límites
de alguna variable, emplear ``None``. Los ``None`` serán remplazados por
el valor (-10**3).
limites_sup : `list` or `numpy.ndarray`
límite superior de cada variable. Si solo se quiere predefinir límites
de alguna variable, emplear ``None``. Los``None`` serán remplazados por
el valor (+10**3).
mejor_valor : `numpy.ndarray`
mejor valor que ha tenido la partícula hasta el momento.
mejor_posicion : `numpy.ndarray`
posición en la que la partícula ha tenido el mejor valor hasta el momento.
valor : `float`
valor actual de la partícula. Resultado de evaluar la función objetivo
con la posición actual.
velocidad : `numpy.ndarray`
array con la velocidad actual de la partícula.
Raises
------
raise Exception
si `limites_inf` es distinto de None y su longitud no coincide con
`n_variables`.
raise Exception
si `limites_sup` es distinto de None y su longitud no coincide con
`n_variables`.
Examples
--------
Ejemplo creación partícula.
>>> part = Particula(
n_variables = 3,
limites_inf = [4,10,20],
limites_sup = [-1,2,0],
verbose = True
)
"""
def __init__(self, n_variables, limites_inf=None, limites_sup=None,
verbose=False):
# Número de variables de la partícula
self.n_variables = n_variables
# Límite inferior de cada variable
self.limites_inf = limites_inf
# Límite superior de cada variable
self.limites_sup = limites_sup
# Posición de la partícula
self.posicion = np.repeat(None, n_variables)
# Velocidad de la parícula
self.velocidad = np.repeat(None, n_variables)
# Valor de la partícula
self.valor = np.repeat(None, 1)
# Mejor valor que ha tenido la partícula hasta el momento
self.mejor_valor = None
# Mejor posición en la que ha estado la partícula hasta el momento
self.mejor_posicion = None
# CONVERSIONES DE TIPO INICIALES
# ----------------------------------------------------------------------
# Si limites_inf o limites_sup no son un array numpy, se convierten en
# ello.
if self.limites_inf is not None \
and not isinstance(self.limites_inf,np.ndarray):
self.limites_inf = np.array(self.limites_inf)
if self.limites_sup is not None \
and not isinstance(self.limites_sup,np.ndarray):
self.limites_sup = np.array(self.limites_sup)
# COMPROBACIONES INICIALES: EXCEPTIONS Y WARNINGS
# ----------------------------------------------------------------------
if self.limites_inf is not None \
and len(self.limites_inf) != self.n_variables:
raise Exception(
"limites_inf debe tener un valor por cada variable. " +
"Si para alguna variable no se quiere límite, emplear None. " +
"Ejemplo: limites_inf = [10, None, 5]"
)
elif self.limites_sup is not None \
and len(self.limites_sup) != self.n_variables:
raise Exception(
"limites_sup debe tener un valor por cada variable. " +
"Si para alguna variable no se quiere límite, emplear None. " +
"Ejemplo: limites_sup = [10, None, 5]"
)
elif (self.limites_inf is None) or (self.limites_sup is None):
warnings.warn(
"Es altamente recomendable indicar los límites dentro de los " +
"cuales debe buscarse la solución de cada variable. " +
"Por defecto se emplea [-10^3, 10^3]."
)
elif any(np.concatenate((self.limites_inf, self.limites_sup)) == None):
warnings.warn(
"Los límites empleados por defecto cuando no se han definido " +
"son: [-10^3, 10^3]."
)
# COMPROBACIONES INICIALES: ACCIONES
# ----------------------------------------------------------------------
# Si no se especifica limites_inf, el valor mínimo que pueden tomar las
# variables es -10^3.
if self.limites_inf is None:
self.limites_inf = np.repeat(-10**3, self.n_variables)
# Si no se especifica limites_sup, el valor máximo que pueden tomar las
# variables es 10^3.
if self.limites_sup is None:
self.limites_sup = np.repeat(+10**3, self.n_variables)
# Si los límites no son nulos, se reemplazan aquellas posiciones None por
# el valor por defecto -10^3 y 10^3.
if self.limites_inf is not None:
self.limites_inf[self.limites_inf == None] = -10**3
if self.limites_sup is not None:
self.limites_sup[self.limites_sup == None] = +10**3
# BUCLE PARA ASIGNAR UN VALOR A CADA UNA DE LAS VARIABLES QUE DEFINEN LA
# POSICIÓN
# ----------------------------------------------------------------------
for i in np.arange(self.n_variables):
# Para cada posición, se genera un valor aleatorio dentro del rango
# permitido para esa variable.
self.posicion[i] = random.uniform(
self.limites_inf[i],
self.limites_sup[i]
)
# LA VELOCIDAD INICIAL DE LA PARTÍCULA ES 0
# ----------------------------------------------------------------------
self.velocidad = np.repeat(0, self.n_variables)
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------
if verbose:
print("Nueva partícula creada")
print("----------------------")
print("Posición: " + str(self.posicion))
print("Límites inferiores de cada variable: " \
+ str(self.limites_inf))
print("Límites superiores de cada variable: " \
+ str(self.limites_sup))
print("Velocidad: " + str(self.velocidad))
print("")
def __repr__(self):
"""
Información que se muestra cuando se imprime un objeto partícula.
"""
texto = "Partícula" \
+ "\n" \
+ "---------" \
+ "\n" \
+ "Posición: " + str(self.posicion) \
+ "\n" \
+ "Velocidad: " + str(self.velocidad) \
+ "\n" \
+ "Mejor posicion: " + str(self.mejor_posicion) \
+ "\n" \
+ "Mejor valor: " + str(self.mejor_valor) \
+ "\n" \
+ "Límites inferiores de cada variable: " \
+ str(self.limites_inf) \
+ "\n" \
+ "Límites superiores de cada variable: " \
+ str(self.limites_sup) \
+ "\n"
return(texto)
def evaluar_particula(self, funcion_objetivo, optimizacion, verbose = False):
"""Este método evalúa una partícula calculando el valor que toma la
función objetivo en la posición en la que se encuentra. Además, compara
si la nueva posición es mejor que las anteriores. Modifica los atributos
valor, mejor_valor y mejor_posicion de la partícula.
Parameters
----------
funcion_objetivo : `function`
función que se quiere optimizar.
optimizacion : {'maximizar', 'minimizar'}
dependiendo de esto, el mejor valor histórico de la partícula será
el mayor o el menor valor que ha tenido hasta el momento.
verbose : `bool`, optional
mostrar información del proceso por pantalla. (default is ``False``)
Raises
------
raise Exception
si el argumento `optimizacion` es distinto de 'maximizar' o 'minimizar'
Examples
--------
Ejemplo evaluar partícula con una función objetivo.
>>> part = Particula(
n_variables = 3,
limites_inf = [4,10,20],
limites_sup = [-1,2,None],
verbose = True
)
>>> def funcion_objetivo(x_0, x_1, x_2):
f= x_0**2 + x_1**2 + x_2**2
return(f)
>>> part.evaluar_particula(
funcion_objetivo = funcion_objetivo,
optimizacion = "maximizar",
verbose = True
)
"""
# COMPROBACIONES INICIALES: EXCEPTIONS Y WARNINGS
# ----------------------------------------------------------------------
if not optimizacion in ["maximizar", "minimizar"]:
raise Exception(
"El argumento optimizacion debe ser: maximizar o minimizar"
)
# EVALUACIÓN DE LA FUNCIÓN OBJETIVO EN LA POSICIÓN ACTUAL
# ----------------------------------------------------------------------
self.valor = funcion_objetivo(*self.posicion)
# MEJOR VALOR Y POSICIÓN
# ----------------------------------------------------------------------
# Se compara el valor actual con el mejor valor histórico. La comparación
# es distinta dependiendo de si se desea maximizar o minimizar.
# Si no existe ningún valor histórico, se almacena el actual. Si ya existe
# algún valor histórico se compara con el actual y, de ser mejor este
# último, se sobrescribe.
if self.mejor_valor is None:
self.mejor_valor = np.copy(self.valor)
self.mejor_posicion = np.copy(self.posicion)
else:
if optimizacion == "minimizar":
if self.valor < self.mejor_valor:
self.mejor_valor = np.copy(self.valor)
self.mejor_posicion = np.copy(self.posicion)
else:
if self.valor > self.mejor_valor:
self.mejor_valor = np.copy(self.valor)
self.mejor_posicion = np.copy(self.posicion)
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------
if verbose:
print("La partícula ha sido evaluada")
print("-----------------------------")
print("Valor actual: " + str(self.valor))
print("")
def mover_particula(self, mejor_p_enjambre, inercia=0.8, peso_cognitivo=2,
peso_social=2, verbose=False):
"""
Este método ejecuta el movimiento de una partícula, lo que implica
actualizar su velocidad y posición. No se permite que la partícula
salga de la zona de búsqueda acotada por los límites.
Parameters
----------
mejor_p_enjambre : `np.narray`
mejor posición de todo el enjambre.
inercia : `float`, optional
coeficiente de inercia. (default is 0.8)
peso_cognitivo : `float`, optional
coeficiente cognitivo. (default is 2)
peso_social : `float`, optional
coeficiente social. (default is 2)
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
Examples
--------
Ejemplo mover partícula.
>>> part = Particula(
n_variables = 3,
limites_inf = [4,10,20],
limites_sup = [-1,2,None],
verbose = True
)
>>> def funcion_objetivo(x_0, x_1, x_2):
f= x_0**2 + x_1**2 + x_2**2
return(f)
>>> part.evaluar_particula(
funcion_objetivo = funcion_objetivo,
optimizacion = "maximizar",
verbose = True
)
>>> part.mover_particula(
mejor_p_enjambre = np.array([-1000,-1000,+1000]),
inercia = 0.8,
peso_cognitivo = 2,
peso_social = 2,
verbose = True
)
"""
# ACTUALIZACIÓN DE LA VELOCIDAD
# ----------------------------------------------------------------------
componente_velocidad = inercia * self.velocidad
r1 = np.random.uniform(low=0.0, high=1.0, size = len(self.velocidad))
r2 = np.random.uniform(low=0.0, high=1.0, size = len(self.velocidad))
componente_cognitivo = peso_cognitivo * r1 * (self.mejor_posicion \
- self.posicion)
componente_social = peso_social * r2 * (mejor_p_enjambre \
- self.posicion)
nueva_velocidad = componente_velocidad + componente_cognitivo \
+ componente_social
self.velocidad = np.copy(nueva_velocidad)
# ACTUALIZACIÓN DE LA POSICIÓN
# ----------------------------------------------------------------------
self.posicion = self.posicion + self.velocidad
# COMPROBAR LÍMITES
# ----------------------------------------------------------------------
# Se comprueba si algún valor de la nueva posición supera los límites
# impuestos. En tal caso, se sobrescribe con el valor del límite
# correspondiente y se reinicia a 0 la velocidad de la partícula en esa
# componente.
for i in np.arange(len(self.posicion)):
if self.posicion[i] < self.limites_inf[i]:
self.posicion[i] = self.limites_inf[i]
self.velocidad[i] = 0
if self.posicion[i] > self.limites_sup[i]:
self.posicion[i] = self.limites_sup[i]
self.velocidad[i] = 0
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------
if verbose:
print("La partícula se ha desplazado")
print("-----------------------------")
print("Nueva posición: " + str(self.posicion))
print("")
Clase enjambre¶
################################################################################
# CLASE ENJAMBRE (SWARM) #
################################################################################
class Enjambre:
"""
Esta clase crea un enjambre de n partículas. El rango de posibles valores
para cada variable (posición) puede estar acotado.
Parameters
----------
n_particulas :`int`
número de partículas del enjambre.
n_variables : `int`
número de variables que definen la posición de las partícula.
limites_inf : `list` or `numpy.ndarray`
límite inferior de cada variable. Si solo se quiere predefinir límites
de alguna variable, emplear ``None``. Los ``None`` serán remplazados por
el valor (-10**3).
limites_sup : `list` or `numpy.ndarray`
límite superior de cada variable. Si solo se quiere predefinir límites
de alguna variable, emplear ``None``. Los``None`` serán remplazados por
el valor (+10**3).
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
Attributes
----------
partículas : `list`
lista con todas las partículas del enjambre.
n_particulas :`int`
número de partículas del enjambre.
n_variables : `int`
número de variables que definen la posición de las partícula.
limites_inf : `list` or `numpy.ndarray`
límite inferior de cada variable.
limites_sup : `list` or `numpy.ndarray`
límite superior de cada variable.
mejor_particula : `object particula`
la mejor partícula del enjambre en estado actual.
mejor_valor : `floar`
el mejor valor del enjambre en su estado actual.
historico_particulas : `list`
lista con el estado de las partículas en cada una de las iteraciones que
ha tenido el enjambre.
historico_mejor_posicion : `list`
lista con la mejor posición en cada una de las iteraciones que ha tenido
el enjambre.
historico_mejor_valor : `list`
lista con el mejor valor en cada una de las iteraciones que ha tenido el
enjambre.
diferencia_abs : `list`
diferencia absoluta entre el mejor valor de iteraciones consecutivas.
resultados_df : `pandas.core.frame.DataFrame`
dataframe con la información del mejor valor y posición encontrado en
cada iteración, así como la mejora respecto a la iteración anterior.
valor_optimo : `float`
mejor valor encontrado en todas las iteraciones.
posicion_optima : `numpy.narray`
posición donde se ha encontrado el valor_optimo.
optimizado : `bool`
si el enjambre ha sido optimizado.
iter_optimizacion : `int`
número de iteraciones de optimizacion.
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
Examples
--------
Ejemplo crear enjambre
>>> enjambre = Enjambre(
n_particulas = 5,
n_variables = 3,
limites_inf = [-5,-5,-5],
limites_sup = [5,5,5],
verbose = True
)
"""
def __init__(self, n_particulas, n_variables, limites_inf = None,
limites_sup = None, verbose = False):
# Número de partículas del enjambre
self.n_particulas = n_particulas
# Número de variables de cada partícula
self.n_variables = n_variables
# Límite inferior de cada variable
self.limites_inf = limites_inf
# Límite superior de cada variable
self.limites_sup = limites_sup
# Verbose
self.verbose = verbose
# Lista de las partículas del enjambre
self.particulas = []
# Etiqueta para saber si el enjambre ha sido optimizado
self.optimizado = False
# Número de iteraciones de optimización llevadas a cabo
self.iter_optimizacion = None
# Mejor partícula del enjambre
self.mejor_particula = None
# Mejor valor del enjambre
self.mejor_valor = None
# Posición del mejor valor del enjambre.
self.mejor_posicion = None
# Estado de todas las partículas del enjambre en cada iteración.
self.historico_particulas = []
# Mejor posición en cada iteración.
self.historico_mejor_posicion = []
# Mejor valor en cada iteración.
self.historico_mejor_valor = []
# Diferencia absoluta entre el mejor valor de iteraciones consecutivas.
self.diferencia_abs = []
# data.frame con la información del mejor valor y posición encontrado en
# cada iteración, así como la mejora respecto a la iteración anterior.
self.resultados_df = None
# Mejor valor de todas las iteraciones
self.valor_optimo = None
# Mejor posición de todas las iteraciones
self.posicion_optima = None
# CONVERSIONES DE TIPO INICIALES
# ----------------------------------------------------------------------
# Si limites_inf o limites_sup no son un array numpy, se convierten en
# ello.
if self.limites_inf is not None \
and not isinstance(self.limites_inf,np.ndarray):
self.limites_inf = np.array(self.limites_inf)
if self.limites_sup is not None \
and not isinstance(self.limites_sup,np.ndarray):
self.limites_sup = np.array(self.limites_sup)
# SE CREAN LAS PARTÍCULAS DEL ENJAMBRE Y SE ALMACENAN
# ----------------------------------------------------------------------
for i in np.arange(n_particulas):
particula_i = Particula(
n_variables = self.n_variables,
limites_inf = self.limites_inf,
limites_sup = self.limites_sup,
verbose = self.verbose
)
self.particulas.append(particula_i)
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------
if verbose:
print("---------------")
print("Enjambre creado")
print("---------------")
print("Número de partículas: " + str(self.n_particulas))
print("Límites inferiores de cada variable: " \
+ np.array2string(self.limites_inf))
print("Límites superiores de cada variable: " \
+ np.array2string(self.limites_sup))
print("")
def __repr__(self):
"""
Información que se muestra cuando se imprime un objeto enjambre.
"""
texto = "============================" \
+ "\n" \
+ " Enjambre" \
+ "\n" \
+ "============================" \
+ "\n" \
+ "Número de partículas: " + str(self.n_particulas) \
+ "\n" \
+ "Límites inferiores de cada variable: " + str(self.limites_inf) \
+ "\n" \
+ "Límites superiores de cada variable: " + str(self.limites_sup) \
+ "\n" \
+ "Optimizado: " + str(self.optimizado) \
+ "\n" \
+ "Iteraciones optimización: " + str(self.iter_optimizacion) \
+ "\n" \
+ "\n" \
+ "Información mejor partícula:" \
+ "\n" \
+ "----------------------------" \
+ "\n" \
+ "Mejor posición actual: " + str(self.mejor_posicion) \
+ "\n" \
+ "Mejor valor actual: " + str(self.mejor_valor) \
+ "\n" \
+ "\n" \
+ "Resultados tras optimizar:" \
+ "\n" \
+ "----------------------------" \
+ "\n" \
+ "Posición óptima: " + str(self.posicion_optima) \
+ "\n" \
+ "Valor óptimo: " + str(self.valor_optimo)
return(texto)
def mostrar_particulas(self, n=None):
"""
Este método muestra la información de cada una de las n primeras
partículas del enjambre.
Parameters
----------
n : `int`
número de particulas que se muestran. Si no se indica el valor
(por defecto ``None``), se muestran todas. Si el valor es mayor
que `self.n_particulas` se muestran todas.
Examples
--------
>>> enjambre = Enjambre(
n_particulas = 5,
n_variables = 3,
limites_inf = [-5,-5,-5],
limites_sup = [5,5,5],
verbose = True
)
>>> enjambre.mostrar_particulas(n = 1)
"""
if n is None:
n = self.n_particulas
elif n > self.n_particulas:
n = self.n_particulas
for i in np.arange(n):
print(self.particulas[i])
return(None)
def evaluar_enjambre(self, funcion_objetivo, optimizacion, verbose = False):
"""
Este método evalúa todas las partículas del enjambre, actualiza sus
valores e identifica la mejor partícula.
Parameters
----------
funcion_objetivo : `function`
función que se quiere optimizar.
optimizacion : {maximizar o minimizar}
Dependiendo de esto, el mejor valor histórico de la partícula será
el mayor o el menor valorque ha tenido hasta el momento.
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
Examples
--------
Ejemplo evaluar enjambre
>>> enjambre = Enjambre(
n_particulas = 5,
n_variables = 3,
limites_inf = [-5,-5,-5],
limites_sup = [5,5,5],
verbose = True
)
>>> def funcion_objetivo(x_0, x_1, x_2):
f= x_0**2 + x_1**2 + x_2**2
return(f)
>>> enjambre.evaluar_enjambre(
funcion_objetivo = funcion_objetivo,
optimizacion = "minimizar",
verbose = True
)
"""
# SE EVALÚA CADA PARTÍCULA DEL ENJAMBRE
# ----------------------------------------------------------------------
for i in np.arange(self.n_particulas):
self.particulas[i].evaluar_particula(
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
verbose = verbose
)
# MEJOR PARTÍCULA DEL ENJAMBRE
# ----------------------------------------------------------------------
# Se identifica la mejor partícula de todo el enjambre. Si se está
# maximizando, la mejor partícula es aquella con mayor valor.
# Lo contrario si se está minimizando.
# Se selecciona inicialmente como mejor partícula la primera.
self.mejor_particula = copy.deepcopy(self.particulas[0])
# Se comparan todas las partículas del enjambre.
for i in np.arange(self.n_particulas):
if optimizacion == "minimizar":
if self.particulas[i].valor < self.mejor_particula.valor:
self.mejor_particula = copy.deepcopy(self.particulas[i])
else:
if self.particulas[i].valor > self.mejor_particula.valor:
self.mejor_particula = copy.deepcopy(self.particulas[i])
# Se extrae la posición y valor de la mejor partícula y se almacenan
# como mejor valor y posición del enjambre.
self.mejor_valor = self.mejor_particula.valor
self.mejor_posicion = self.mejor_particula.posicion
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------
if verbose:
print("-----------------")
print("Enjambre evaluado")
print("-----------------")
print("Mejor posición encontrada : "
+ np.array2string(self.mejor_posicion))
print("Mejor valor encontrado : " + str(self.mejor_valor))
print("")
def mover_enjambre(self, inercia, peso_cognitivo, peso_social,
verbose = False):
"""
Este método mueve todas las partículas del enjambre.
Parameters
----------
optimizacion : {maximizar o minimizar}
si se desea maximizar o minimizar la función.
inercia : `float` or `int`
coeficiente de inercia.
peso_cognitivo : `float` or `int`
coeficiente cognitivo.
peso_social : `float` or `int`
coeficiente social.
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
"""
# Se actualiza la posición de cada una de las partículas que forman el
# enjambre.
for i in np.arange(self.n_particulas):
self.particulas[i].mover_particula(
mejor_p_enjambre = self.mejor_posicion,
inercia = inercia,
peso_cognitivo = peso_cognitivo,
peso_social = peso_social,
verbose = verbose
)
# Información del proceso (VERBOSE)
# ----------------------------------------------------------------------
if verbose:
print("---------------------------------------------------------" \
"------------")
print("La posición de todas las partículas del enjambre ha sido " \
"actualizada.")
print("---------------------------------------------------------" \
"------------")
print("")
def optimizar(self, funcion_objetivo, optimizacion, n_iteraciones = 50,
inercia = 0.8, reduc_inercia = True, inercia_max = 0.9,
inercia_min = 0.4, peso_cognitivo = 2, peso_social = 2,
parada_temprana = False, rondas_parada = None,
tolerancia_parada = None, verbose = False):
"""
Este método realiza el proceso de optimización de un enjambre.
Parameters
----------
funcion_objetivo : `function`
función que se quiere optimizar.
optimizacion : {maximizar o minimizar}
si se desea maximizar o minimizar la función.
m_iteraciones : `int` , optional
numero de iteraciones de optimización. (default is ``50``)
inercia : `float` or `int`, optional
coeficiente de inercia. (default is ``0.8``)
peso_cognitivo : `float` or `int`, optional
coeficiente cognitivo. (default is ``2``)
peso_social : `float` or `int`, optional
coeficiente social. (default is ``2``)
reduc_inercia: `bool`, optional
activar la reducción del coeficiente de inercia. En tal caso, el
argumento `inercia` es ignorado. (default is ``True``)
inercia_max : `float` or `int`, optional
valor inicial del coeficiente de inercia si se activa `reduc_inercia`.
(default is ``0.9``)
inercia_min : `float` or `int`, optional
valor minimo del coeficiente de inercia si se activa `reduc_min`.
(default is ``0.4``)
parada_temprana : `bool`, optional
si durante las últimas `rondas_parada` generaciones la diferencia
absoluta entre mejores individuos no es superior al valor de
`tolerancia_parada`, se detiene el algoritmo y no se crean nuevas
generaciones. (default is ``False``)
rondas_parada : `int`, optional
número de generaciones consecutivas sin mejora mínima para que se
active la parada temprana. (default is ``None``)
tolerancia_parada : `float` or `int`, optional
valor mínimo que debe tener la diferencia de generaciones consecutivas
para considerar que hay cambio. (default is ``None``)
verbose : `bool`, optional
mostrar información de la partícula creada. (default is ``False``)
Raises
------
raise Exception
si se indica `parada_temprana = True` y los argumentos `rondas_parada`
o `tolerancia_parada` son ``None``.
raise Exception
si se indica `reduc_inercia = True` y los argumentos `inercia_max`
o `inercia_min` son ``None``.
Examples
--------
Ejemplo optimización
>>> def funcion_objetivo(x_0, x_1):
# Para la región acotada entre −10<=x_0<=0 y −6.5<=x_1<=0 la
# función tiene múltiples mínimos locales y un único minimo
# global en f(−3.1302468,−1.5821422)= −106.7645367.
f = np.sin(x_1)*np.exp(1-np.cos(x_0))**2 \
+ np.cos(x_0)*np.exp(1-np.sin(x_1))**2 \
+ (x_0-x_1)**2
return(f)
>>> enjambre = Enjambre(
n_particulas = 50,
n_variables = 2,
limites_inf = [-10, -6.5],
limites_sup = [0, 0],
verbose = False
)
>>> enjambre.optimizar(
funcion_objetivo = funcion_objetivo,
optimizacion = "minimizar",
n_iteraciones = 250,
inercia = 0.8,
reduc_inercia = True,
inercia_max = 0.9,
inercia_min = 0.4,
peso_cognitivo = 1,
peso_social = 2,
parada_temprana = True,
rondas_parada = 5,
tolerancia_parada = 10**-3,
verbose = False
)
"""
# COMPROBACIONES INICIALES: EXCEPTIONS Y WARNINGS
# ----------------------------------------------------------------------
# Si se activa la parada temprana, hay que especificar los argumentos
# rondas_parada y tolerancia_parada.
if parada_temprana \
and (rondas_parada is None or tolerancia_parada is None):
raise Exception(
"Para activar la parada temprana es necesario indicar un " \
+ " valor de rondas_parada y de tolerancia_parada."
)
# Si se activa la reducción de inercia, hay que especificar los argumentos
# inercia_max y inercia_min.
if reduc_inercia \
and (inercia_max is None or inercia_min is None):
raise Exception(
"Para activar la reducción de inercia es necesario indicar un " \
+ "valor de inercia_max y de inercia_min."
)
# ITERACIONES
# ----------------------------------------------------------------------
start = time.time()
for i in np.arange(n_iteraciones):
if verbose:
print("-------------")
print("Iteracion: " + str(i))
print("-------------")
# EVALUAR PARTÍCULAS DEL ENJAMBRE
# ------------------------------------------------------------------
enjambre.evaluar_enjambre(
funcion_objetivo = funcion_objetivo,
optimizacion = "minimizar",
verbose = verbose
)
# SE ALMACENA LA INFORMACIÓN DE LA ITERACIÓN EN LOS HISTÓRICOS
# ------------------------------------------------------------------
self.historico_particulas.append(copy.deepcopy(self.particulas))
self.historico_mejor_posicion.append(copy.deepcopy(self.mejor_posicion))
self.historico_mejor_valor.append(copy.deepcopy(self.mejor_valor))
# SE CALCULA LA DIFERENCIA ABSOLUTA RESPECTO A LA ITERACION ANTERIOR
# ------------------------------------------------------------------
# La diferencia solo puede calcularse a partir de la segunda
# iteración.
if i == 0:
self.diferencia_abs.append(None)
else:
diferencia = abs(self.historico_mejor_valor[i] \
- self.historico_mejor_valor[i-1])
self.diferencia_abs.append(diferencia)
# CRITERIO DE PARADA
# ------------------------------------------------------------------
# Si durante las últimas n generaciones, la diferencia absoluta entre
# mejores partículas no es superior al valor de tolerancia_parada,
# se detiene el algoritmo y no se crean nuevas generaciones.
if parada_temprana and i > rondas_parada:
ultimos_n = np.array(self.diferencia_abs[-(rondas_parada): ])
if all(ultimos_n < tolerancia_parada):
print("Algoritmo detenido en la iteracion "
+ str(i) \
+ " por falta cambio absoluto mínimo de " \
+ str(tolerancia_parada) \
+ " durante " \
+ str(rondas_parada) \
+ " iteraciones consecutivas.")
break
# MOVER PARTÍCULAS DEL ENJAMBRE
# ------------------------------------------------------------------
# Si se ha activado la reducción de inercia, se recalcula su valor
# para la iteración actual.
if reduc_inercia:
inercia = ((inercia_max - inercia_min) \
* (n_iteraciones-i)/n_iteraciones) \
+ inercia_min
enjambre.mover_enjambre(
inercia = inercia,
peso_cognitivo = peso_cognitivo,
peso_social = peso_social,
verbose = False
)
end = time.time()
self.optimizado = True
self.iter_optimizacion = i
# IDENTIFICACIÓN DEL MEJOR INDIVIDUO DE TODO EL PROCESO
# ----------------------------------------------------------------------
indice_valor_optimo = np.argmin(np.array(self.historico_mejor_valor))
self.valor_optimo = self.historico_mejor_valor[indice_valor_optimo]
self.posicion_optima = self.historico_mejor_posicion[indice_valor_optimo]
# CREACIÓN DE UN DATAFRAME CON LOS RESULTADOS
# ----------------------------------------------------------------------
self.resultados_df = pd.DataFrame(
{
"mejor_valor_enjambre" : self.historico_mejor_valor,
"mejor_posicion_enjambre": self.historico_mejor_posicion,
"diferencia_abs" : self.diferencia_abs
}
)
self.resultados_df["iteracion"] = self.resultados_df.index
print("-------------------------------------------")
print("Optimización finalizada " \
+ datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
print("-------------------------------------------")
print("Duración optimización: " + str(end - start))
print("Número de iteraciones: " + str(self.iter_optimizacion))
Ejemplos¶
Crear partícula¶
# Ejemplo creación partícula.
part = Particula(
n_variables = 3,
limites_inf = [4,10,20],
limites_sup = [-1,2,5],
verbose = True
)
part
Evaluar partícula¶
# Ejemplo evaluar partícula con una función objetivo
def funcion_objetivo(x_0, x_1, x_2):
f= x_0**2 + x_1**2 + x_2**2
return(f)
part.evaluar_particula(
funcion_objetivo = funcion_objetivo,
optimizacion = "maximizar",
verbose = True
)
# Hasta que la partícula se mueva, el valor actual y mejor valor es el mismo.
part
Mover partícula¶
# Ejemplo mover partícula
part.mover_particula(
mejor_p_enjambre = np.array([-1000,-1000,+1000]),
inercia = 0.8,
peso_cognitivo = 2,
peso_social = 2,
verbose = True
)
part
Crear enjambre¶
# Ejemplo crear enjambre
enjambre = Enjambre(
n_particulas = 4,
n_variables = 3,
limites_inf = [-5,-5,-5],
limites_sup = [5,5,5],
verbose = True
)
Evaluar enjambre¶
# Ejemplo evaluar enjambre
def funcion_objetivo(x_0, x_1, x_2):
f= x_0**2 + x_1**2 + x_2**2
return(f)
enjambre.evaluar_enjambre(
funcion_objetivo = funcion_objetivo,
optimizacion = "minimizar",
verbose = True
)
Mover enjambre¶
# Ejemplo mover enjambre
enjambre.mover_enjambre(
inercia = 0.8,
peso_cognitivo = 2,
peso_social = 2,
verbose = True
)
Optimización¶
# Ejemplo optimización
def funcion_objetivo(x_0, x_1):
'''
Para la región acotada entre −10<=x_0<=0 y −6.5<=x_1<=0 la función tiene
múltiples mínimos locales y un único minimo global que se encuentra en
f(−3.1302468,−1.5821422) = −106.7645367
'''
f = np.sin(x_1)*np.exp(1-np.cos(x_0))**2 \
+ np.cos(x_0)*np.exp(1-np.sin(x_1))**2 \
+ (x_0-x_1)**2
return(f)
# Contour plot función objetivo
x_0 = np.linspace(start = -10, stop = 0, num = 100)
x_1 = np.linspace(start = -6.5, stop = 0, num = 100)
x_0, x_1 = np.meshgrid(x_0, x_1)
z = funcion_objetivo(x_0, x_1)
plt.contour(x_0, x_1, z, 35, cmap='RdGy');
enjambre = Enjambre(
n_particulas = 50,
n_variables = 2,
limites_inf = [-10, -6.5],
limites_sup = [0, 0],
verbose = False
)
enjambre.optimizar(
funcion_objetivo = funcion_objetivo,
optimizacion = "minimizar",
n_iteraciones = 250,
inercia = 0.8,
reduc_inercia = True,
inercia_max = 0.9,
inercia_min = 0.4,
peso_cognitivo = 1,
peso_social = 2,
parada_temprana = True,
rondas_parada = 5,
tolerancia_parada = 10**-3,
verbose = False
)
enjambre
# Evolución de la optimización
fig = plt.figure(figsize=(6,4))
enjambre.resultados_df['mejor_valor_enjambre'].plot();
# Representación evolución partículas gráfico animado.
# Se extrae la posición de las partículas en cada iteración del enjambre
import plotly_express as px
def extraer_posicion(particula):
posicion = particula.posicion
return(posicion)
lista_df_temp = []
for i in np.arange(len(enjambre.historico_particulas)):
posiciones = list(map(extraer_posicion, enjambre.historico_particulas[i]))
df_temp = pd.DataFrame({"iteracion": i, "posicion": posiciones})
lista_df_temp.append(df_temp)
df_posiciones = pd.concat(lista_df_temp)
df_posiciones[['x_0','x_1']] = pd.DataFrame(df_posiciones["posicion"].values.tolist(),
index= df_posiciones.index)
df_posiciones.head()
px.scatter(
df_posiciones,
x = "x_0",
y = "x_1",
range_x = [-10, 0],
range_y = [-6.5, 0],
animation_frame = "iteracion"
)
#import matplotlib.animation as animation
#%matplotlib notebook
#fig = plt.figure(figsize=(8,5))
#plt.xlim(-10,0)
#plt.ylim(-6.5,0)
#def animate(i):
# p2 = fig.clear()
# plt.xlim(-10,0)
# plt.ylim(-6.5,0)
# df_posiciones_i = df_posiciones[df_posiciones["iteracion"] == i][["x_0", "x_1"]] #select data range
# p1 = plt.contour(x_0, x_1, z, 35, cmap='RdGy')
# p2 = plt.scatter(df_posiciones_i["x_0"], df_posiciones_i["x_1"])
#ani = matplotlib.animation.FuncAnimation(fig, animate, repeat = True, blit = True)
# Guardar la animación como mp4
#Writer = animation.writers['ffmpeg']
#writer = Writer(fps=8, bitrate=1800)
#ani.save('animacion_pso.mp4', writer=writer)
%%HTML
<div align="middle">
<video width="60%" controls>
<source src="./animacion_pso.mp4" type="video/mp4">
</video></div>
Información de sesión¶
from sinfo import sinfo
sinfo()
Bibliografía¶
Particle Swarm Optimization: A TutorialJames BlondinSeptember 4, 2009
Sengupta, Saptarshi & Basak, Sanchita & Alan Peters II, Richard. (2018). Particle Swarm Optimization: A survey of historical and recent developments with hybridization perspectives. 10.3390/make1010010.
Bonyadi, Mohammad reza & Michalewicz, Zbigniew & Li, Xiaodong. (2014). An analysis of the velocity updating rule of the particle swarm optimization algorithm. Journal of Heuristics.
Particle Swarm Optimization: A TutorialJames BlondinSeptember 4, 2009
wikipedia.org/wiki/Particle_swarm_optimization.
¿Cómo citar este documento?
Optimización con enjambre de partículas (Particle Swarm Optimization) por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py02_optimizacion_pso
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
-
Compartir: copiar y redistribuir el material en cualquier medio o formato.
-
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
-
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.