En primer lugar, definamos nuestros hipoparámetros. Como en muchos otros algoritmos metaheurísticos, estas variables deben ajustarse sobre la marcha y no existe un conjunto de valores versátil. Pero ciñámonos a estos:
POP_SIZE = 10 #population size
MAX_ITER = 30 #the amount of optimization iterations
w = 0.2 #inertia weight
c1 = 1 #personal acceleration factor
c2 = 2 #social acceleration factor
Ahora creemos una función que generaría una población aleatoria:
def populate(size):
x1,x2 = -10, 3 #x1, x2 = right and left boundaries of our X axis
pop = rnd.uniform(x1,x2, size) # size = amount of particles in population
return pop
Si lo visualizamos, obtendremos algo como esto:
x1=populate(50)
y1=function(x1)plt.plot(x,y, lw=3, label='Func to optimize')
plt.plot(x1,y1,marker='o', ls='', label='Particles')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()
Aquí puedes ver que inicialicé aleatoriamente una población de 50 partículas, algunas de las cuales ya están cerca de la solución.
Ahora implementemos el propio algoritmo PSO. Comenté cada fila del código, pero si tiene alguna pregunta, no dude en hacerla en los comentarios a continuación.
"""Particle Swarm Optimization (PSO)"""
particles = populate(POP_SIZE) #generating a set of particles
velocities = np.zeros(np.shape(particles)) #velocities of the particles
gains = -np.array(function(particles)) #calculating function values for the populationbest_positions = np.copy(particles) #it's our first iteration, so all positions are the best
swarm_best_position = particles[np.argmax(gains)] #x with with the highest gain
swarm_best_gain = np.max(gains) #highest gain
l = np.empty((MAX_ITER, POP_SIZE)) #array to collect all pops to visualize afterwards
for i in range(MAX_ITER):
l[i] = np.array(np.copy(particles)) #collecting a pop to visualize
r1 = rnd.uniform(0, 1, POP_SIZE) #defining a random coefficient for personal behavior
r2 = rnd.uniform(0, 1, POP_SIZE) #defining a random coefficient for social behavior
velocities = np.array(w * velocities + c1 * r1 * (best_positions - particles) + c2 * r2 * (swarm_best_position - particles)) #calculating velocities
particles+=velocities #updating position by adding the velocity
new_gains = -np.array(function(particles)) #calculating new gains
idx = np.where(new_gains > gains) #getting index of Xs, which have a greater gain now
best_positions[idx] = particles[idx] #updating the best positions with the new particles
gains[idx] = new_gains[idx] #updating gains
if np.max(new_gains) > swarm_best_gain: #if current maxima is greateer than across all previous iters, than assign
swarm_best_position = particles[np.argmax(new_gains)] #assigning the best candidate solution
swarm_best_gain = np.max(new_gains) #assigning the best gain
print(f'Iteration {i+1} \tGain: {swarm_best_gain}')
Después de 30 iteraciones tenemos esto:
Como puede ver, el algoritmo cayó en el mínimo local, que no es lo que queríamos. Por eso necesitamos ajustar nuestros hipoparámetros y empezar de nuevo. Esta vez decidí poner peso de inercia. p=0,8por lo tanto, ahora la velocidad anterior tiene un mayor impacto en el estado actual.
Y listo, llegamos al mínimo global de la función. Te recomiendo encarecidamente que juegues con POP_SIZE, c₁ y c₂. Le permitirá comprender mejor el código y la idea detrás de PSO. Si está interesado, puede complicar la tarea y optimizar alguna función 3D y hacer una buena visualización.
=============================================
[1]Shi Y. Optimización de enjambre de partículas // Conexiones IEEE. — 2004. — Т. 2. — №. 1. — С. 8–13.
=============================================
Todos mis artículos en Medium son gratuitos y de acceso abierto, ¡por eso te agradecería mucho que me siguieras aquí!
PD: Soy un apasionado de la ciencia (geo)datos, el aprendizaje automático/inteligencia artificial y el cambio climático. Entonces, si quieres trabajar juntos en algún proyecto, por favor contáctame en LinkedIn.
🛰️Sigue para más🛰️