0ufjlaenxgswppj8s.jpeg

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()

Imagen por autor.

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 population

best_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:

PSO (w=0,2, c1=1, c2=2). Imagen por autor.

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.

PSO (w=0,9, c1=1, c2=2). Imagen por autor.

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🛰️