# Particle Swarm Optimization Algorithm (explained with Python!)

[SPOILER] We will be using the Particle Swarm Optimization algorithm to obtain the minumum of some test functions

First of all, let's import the libraries we'll need (remember we are using Python 3)

``````

In [1]:

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from plotPSO import plotPSO_2D
import optitestfuns

# import scipy as sp
# import time

%matplotlib inline
plt.style.use('bmh')

``````

We can define and plot the function we want to optimize:

``````

In [2]:

# Testing 2D plot

lo_b = -5 # lower bound
up_b =  5 # upper bound

limits=([lo_b, up_b], # x bounds
[lo_b, up_b]) # y bounds

x_lo = limits[0][0]
x_up = limits[0][1]
y_lo = limits[1][0]
y_up = limits[1][1]

f = optitestfuns.ackley # Objective function (aka fitness function)

#fig, ax = plotPSO_2D(f,  limits)

``````

## PSO Algorithm

``````

In [3]:

n_iterations = 50

def run_PSO(n_particles=10, omega=0.3, phi_p=0.7, phi_g=0.7):
""" PSO algorithm to a funcion already defined.
Params:
omega = 0.3  # Particle weight (intertial)
phi_p = 0.7  # particle best weight
phi_g = 0.7  # global global weight

"""
global x_best_p_global, y_best_p_global, z_p_best_global, \
x_particles, y_particles, z_particles, \
u_particles, v_particles

# Note: we are using global variables to ease the use of interactive widgets
# This code will work fine without the global (and actually it will be safer)

# Initialazing x postion of particles

x_particles = np.zeros((n_iterations, n_particles))
x_particles[0,:] = np.random.uniform(lo_b, up_b, size=n_particles)

# Initialazing y postion of particles
y_particles = np.zeros((n_iterations, n_particles))
y_particles[0,:] = np.random.uniform(lo_b, up_b, size=n_particles)

# Initialazing best praticles
x_best_particles = np.copy(x_particles[0,:])
y_best_particles = np.copy(y_particles[0,:])

# Calculate Objective function (aka fitness function)
z_particles = np.zeros((n_iterations, n_particles))

for i in range(n_particles):
z_particles[0,i] = f((x_particles[0,i],y_particles[0,i]))

z_best_global = np.min(z_particles[0,:])
index_best_global = np.argmin(z_particles[0,:])

x_best_p_global = x_particles[0,index_best_global]
y_best_p_global = y_particles[0,index_best_global]

# Initialazin velocity
velocity_lo = lo_b-up_b  # [L/iteration]
velocity_up = up_b-lo_b  # [L/iteration]

v_max = 0.07 # [L/iteration]

u_particles = np.zeros((n_iterations, n_particles))
u_particles[0,:] = 0.1*np.random.uniform(velocity_lo, velocity_up, size=n_particles)

v_particles = np.zeros((n_iterations, n_particles))
v_particles[0,:] = 0.1*np.random.uniform(velocity_lo, velocity_up, size=n_particles)

# PSO STARTS
iteration = 1
while iteration <= n_iterations-1:
for i in range(n_particles):
x_p = x_particles[iteration-1, i]
y_p = y_particles[iteration-1, i]

u_p = u_particles[iteration-1, i]
v_p = v_particles[iteration-1, i]

x_best_p = x_best_particles[i]
y_best_p = y_best_particles[i]

r_p = np.random.uniform(0, 1)
r_g = np.random.uniform(0, 1)

u_p_new = omega*u_p + \
phi_p*r_p*(x_best_p-x_p) + \
phi_g*r_g*(x_best_p_global-x_p)

v_p_new = omega*v_p + \
phi_p*r_p*(y_best_p-y_p) + \
phi_g*r_g*(y_best_p_global-y_p)

# # Velocity control
# while not (-v_max <= u_p_new <= v_max):
#     u_p_new = 0.9*u_p_new
# while not (-v_max <= u_p_new <= v_max):
#     u_p_new = 0.9*u_p_new

x_p_new = x_p + u_p_new
y_p_new = y_p + v_p_new

# Ignore new position if it's out of the domain
if not ((lo_b <= x_p_new <= up_b) and (lo_b <= y_p_new <= up_b)):
x_p_new = x_p
y_p_new = y_p

x_particles[iteration, i] = x_p_new
y_particles[iteration, i] = y_p_new

u_particles[iteration, i] = u_p_new
v_particles[iteration, i] = v_p_new

# Evaluation
z_p_new = f((x_p_new,  y_p_new))
z_p_best = f((x_best_p, y_best_p))

z_particles[iteration, i] = z_p_new

if z_p_new < z_p_best:
x_best_particles[i] = x_p_new
y_best_particles[i] = y_p_new

z_p_best_global = f([x_best_p_global, y_best_p_global])

if z_p_new < z_p_best_global:
x_best_p_global = x_p_new
y_best_p_global = y_p_new

# end while loop particles
iteration = iteration + 1

# Plotting convergence
z_particles_best_hist = np.min(z_particles, axis=1)
z_particles_worst_hist = np.max(z_particles, axis=1)

z_best_global = np.min(z_particles)
index_best_global = np.argmin(z_particles)

fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(10, 2))

# Grid points
x_lo = limits[0][0]
x_up = limits[0][1]
y_lo = limits[1][0]
y_up = limits[1][1]

assert x_lo<x_up, "Unbound x limits, the first value of the list needs to be higher"
assert y_lo<y_up, "Unbound x limits, the first value of the list needs to be higher"

n_points = 100

x = np.linspace(x_lo, x_up, n_points) # x coordinates of the grid
y = np.linspace(y_lo, y_up, n_points) # y coordinates of the grid

XX, YY = np.meshgrid(x,y)
ZZ = np.zeros_like(XX)

for i in range(n_points):
for j in range(n_points):
ZZ[i,j] = f((XX[i,j], YY[i,j]))

# Limits of the function being plotted
ax1.plot((0,n_iterations),(np.min(ZZ),np.min(ZZ)), '--g', label="min\$f(x)\$")
ax1.plot((0,n_iterations),(np.max(ZZ),np.max(ZZ)),'--r',  label="max\$f(x)\$")

# Convergence of the best particle and worst particle value
ax1.plot(np.arange(n_iterations),z_particles_best_hist,'b',  label="\$p_{best}\$")
ax1.plot(np.arange(n_iterations),z_particles_worst_hist,'k', label="\$p_{worst}\$")

ax1.set_xlim((0,n_iterations))

ax1.set_ylabel('\$f(x)\$')
ax1.set_xlabel('\$i\$ (iteration)')
ax1.set_title('Convergence')

ax1.legend()

``````
``````

In [4]:

run_PSO()

``````
``````

``````

# Animation

``````

In [5]:

from __future__ import print_function
import ipywidgets as widgets
from IPython.display import display, HTML

``````
``````

In [6]:

def plotPSO_iter(i=0): #iteration
"""Visualization of particles and obj. function"""

fig, (ax1, ax2) = plotPSO_2D(f, limits,
particles_xy=(x_particles[i, :],y_particles[i, :]),
particles_uv=(u_particles[i, :],v_particles[i, :]))

``````
``````

In [7]:

w_arg_PSO = widgets.interact_manual(run_PSO,
n_particles=(2,50),
omega=(0,1,0.001),
phi_p=(0,1,0.001),
phi_g=(0,1,0.001),
continuous_update=False)

``````
``````

``````
``````

In [8]:

w_viz_PSO = widgets.interact_manual(plotPSO_iter, i=(0,n_iterations-1), continuous_update=False)

``````
``````

``````

Let's have a look to some examples in case you can't play with the sliders: