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