# Particle Swarm Optimization (PSO) from scratch. Simplest explanation in python

#### How to implement PSO

Before talking about *swarms* and *particles*, let’s briefly discuss **optimization** itself. Basically, optimization is the process of finding the minima or maxima of some function. For instance, when you need to get to your office ASAP and think about which way is the fastest, you’re optimizing your route (in this case it’s a function). In math, there are literally hundreds of ways of optimization, and among them a sub-group called nature-inspired exists.

**Nature-inspired algorithms** are based on phenomena which draw inspiration from natural phenomena or processes. Among the most popular ones are Genetic Algorithm, Cuckoo Search, Ant Colony and **Particle Swarm Optimization ****[1]**** or PSO**.

This tutorial is implemented in python using only *numpy* and *matplotlib*. To follow up you can use this notebook.

### Let’s start with creating a function which we’ll be optimizing using PSO.

If someone asks me to think about functions, the first thing that comes to my mind (it almost reflexive😂) is parabola. So let’s plot it firstly.

import numpy as np

import numpy.random as rnd

import matplotlib.pyplot as plt

x = np.arange(-10,10, 0.01)

y = x**2

plt.plot(x,y, lw=3)

plt.xlabel('x')

plt.ylabel('y')

plt.grid(True)

plt.show()

In case of this function, y=x², the goal of optimization is to find the point x=0, so y has the lowest value possible — 0. But it’s too easy. Let’s overcomplicate everything, but stay in **2D** space. I decided to use the following system of equations:

Or in code:

def function(x):

y = []

for i in x:

if i>-3:

y.append(i**2-4*i+i**3)

else:

y.append(0.2*i**2)

return y

x = np.arange(-10,3, 0.001)

y = function(x)

plt.plot(x,y, lw=3)

plt.title('Function to optimize')

plt.xlabel('x')

plt.ylabel('y')

plt.grid(True)

plt.show()

As you can see I intentionally made two minima: global (left) and local (right). Often when optimizing a function, we need to find the lowest (greatest) value possible, i.e. global minimum (maximum), so I wanted the PSO algorithm to face the challenge of two extrema.

### Now let’s discuss the algorithm itself.

The underlying idea is based on two terms: **population** and **swarm**. Swarm consists of populations which comprise particles.

### particle — — — > population — — → swarm

What is a particle in this case then? A particle is a possible solution. Imagine that you’re investigating a crime and a bunch of suspects are standing in front of you. All of them are “candidates” to be convicted. But in the case of particles, the candidates are cooperating and have social influence.

The particles themselves has two features — **position** and **velocity **(or speed; our particles are flying, right?). So these particles iteratively, one step at a time, change their positions. Their velocity is defined by the following formula:

As you can see, there are three components. The first is **inertia**. This is the velocity of the particle from the previous iteration multiplied by an arbitrary weight **w**. So it basically defines how much the previous speed of the particle affects the current one.

The second one is **personal influence**. It has its own coefficients: an arbitrary weight *c*₁ and random number *r*₁. In the brackets we have a difference between the **best position** of a particle in a **population **and the position of the current particle (**p**).

The third part is **social influence**. It has similar arbitrary weight *c*₂ and random number *r*₂, as well as the difference between the **best position **of a particle in a **swarm** and the position of the current particle (**p**).

After calculating the velocity we update the current positions of the particle simply by summation:

By now there is one big question left:what are these best position in the swarm and population?

To calculate them we need a gain/reward function indicating which solution is closer to the minimum (maximum). In our example, this gain function is the function we plotted. So the particles are Xs, and by substituting each X to the function, we can figure out which X gives the lowest (greatest) value of Y.

Thus, the best position in a population is X, which gives the lowest (greatest) Y on **the current iteration**. And the best position in a swarm is X which gives the lowest (greatest) Y across all **the previous iterations.**

*So simplifying everything, we can say that the idea of the algorithm is the following:*

A bunch of particles, each having a certain position and velocity, are flying together in a search of the global minimum (maximum), forming a population. Populations appear iteratively and live during only one iteration, but they exchange information with each other, so each following population is closer to finding the solution than the previous one.

**Now let’s get back to coding and implement in python what we’ve discussed.**

First of all, let’s define our hypoparameters. Like in many other metaheuristic algorithms, these variables should be adjusted on the way, and there is no versatile set of values. But let’s stick to these ones:

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

Now let’s create a function which would generate a random population:

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

If we visualize it, we’ll get something like this:

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

Here you can see that I randomly initialized a population of 50 particles, some of which are already close to the solution.

Now let’s implement the PSO algorithm itself. I commented each row in the code, but if you have any questions, feel free to ask in the comments below.

"""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}')

After 30 iteration we’ve got this:

As you can see the algorithm fell into the local minimum, which is not what we wanted. That’s why we need to tune our hypoparameters and start again. This time I decided to set inertia weight **w=0.8**, thus, now the previous velocity has a greater impact on the current state.

And voila, we reached the global minimum of the function. I strongly encourage you to play around with POP_SIZE, c₁ and c₂. It’ll allow you to gain a better understanding of the code and the idea behind PSO. If you’re interested you can complicate the task and optimize some 3D function and make a nice visualization.

===========================================

[1]Shi Y. Particle swarm optimization //IEEE connections. — 2004. — Т. 2. — №. 1. — С. 8–13.

===========================================

*All my articles on Medium are free and open-access, that’s why I’d really appreciate if you followed me here!*

P.s. I’m extremely passionate about (Geo)Data Science, ML/AI and Climate Change. So if you want to work together on some project pls contact me in LinkedIn.

🛰️Follow for more🛰️

Particle Swarm Optimization (PSO) from scratch. Simplest explanation in python was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.

#### How to implement PSO

Before talking about *swarms* and *particles*, let’s briefly discuss **optimization** itself. Basically, optimization is the process of finding the minima or maxima of some function. For instance, when you need to get to your office ASAP and think about which way is the fastest, you’re optimizing your route (in this case it’s a function). In math, there are literally hundreds of ways of optimization, and among them a sub-group called nature-inspired exists.

**Nature-inspired algorithms** are based on phenomena which draw inspiration from natural phenomena or processes. Among the most popular ones are Genetic Algorithm, Cuckoo Search, Ant Colony and **Particle Swarm Optimization ****[1]**** or PSO**.

This tutorial is implemented in python using only *numpy* and *matplotlib*. To follow up you can use this notebook.

### Let’s start with creating a function which we’ll be optimizing using PSO.

If someone asks me to think about functions, the first thing that comes to my mind (it almost reflexive😂) is parabola. So let’s plot it firstly.

import numpy as np

import numpy.random as rnd

import matplotlib.pyplot as plt

x = np.arange(-10,10, 0.01)

y = x**2

plt.plot(x,y, lw=3)

plt.xlabel('x')

plt.ylabel('y')

plt.grid(True)

plt.show()

In case of this function, y=x², the goal of optimization is to find the point x=0, so y has the lowest value possible — 0. But it’s too easy. Let’s overcomplicate everything, but stay in **2D** space. I decided to use the following system of equations:

Or in code:

def function(x):

y = []

for i in x:

if i>-3:

y.append(i**2-4*i+i**3)

else:

y.append(0.2*i**2)

return y

x = np.arange(-10,3, 0.001)

y = function(x)

plt.plot(x,y, lw=3)

plt.title('Function to optimize')

plt.xlabel('x')

plt.ylabel('y')

plt.grid(True)

plt.show()

As you can see I intentionally made two minima: global (left) and local (right). Often when optimizing a function, we need to find the lowest (greatest) value possible, i.e. global minimum (maximum), so I wanted the PSO algorithm to face the challenge of two extrema.

### Now let’s discuss the algorithm itself.

The underlying idea is based on two terms: **population** and **swarm**. Swarm consists of populations which comprise particles.

### particle — — — > population — — → swarm

What is a particle in this case then? A particle is a possible solution. Imagine that you’re investigating a crime and a bunch of suspects are standing in front of you. All of them are “candidates” to be convicted. But in the case of particles, the candidates are cooperating and have social influence.

The particles themselves has two features — **position** and **velocity **(or speed; our particles are flying, right?). So these particles iteratively, one step at a time, change their positions. Their velocity is defined by the following formula:

As you can see, there are three components. The first is **inertia**. This is the velocity of the particle from the previous iteration multiplied by an arbitrary weight **w**. So it basically defines how much the previous speed of the particle affects the current one.

The second one is **personal influence**. It has its own coefficients: an arbitrary weight *c*₁ and random number *r*₁. In the brackets we have a difference between the **best position** of a particle in a **population **and the position of the current particle (**p**).

The third part is **social influence**. It has similar arbitrary weight *c*₂ and random number *r*₂, as well as the difference between the **best position **of a particle in a **swarm** and the position of the current particle (**p**).

After calculating the velocity we update the current positions of the particle simply by summation:

By now there is one big question left:what are these best position in the swarm and population?

To calculate them we need a gain/reward function indicating which solution is closer to the minimum (maximum). In our example, this gain function is the function we plotted. So the particles are Xs, and by substituting each X to the function, we can figure out which X gives the lowest (greatest) value of Y.

Thus, the best position in a population is X, which gives the lowest (greatest) Y on **the current iteration**. And the best position in a swarm is X which gives the lowest (greatest) Y across all **the previous iterations.**

*So simplifying everything, we can say that the idea of the algorithm is the following:*

A bunch of particles, each having a certain position and velocity, are flying together in a search of the global minimum (maximum), forming a population. Populations appear iteratively and live during only one iteration, but they exchange information with each other, so each following population is closer to finding the solution than the previous one.

**Now let’s get back to coding and implement in python what we’ve discussed.**

First of all, let’s define our hypoparameters. Like in many other metaheuristic algorithms, these variables should be adjusted on the way, and there is no versatile set of values. But let’s stick to these ones:

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

Now let’s create a function which would generate a random population:

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

If we visualize it, we’ll get something like this:

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

Here you can see that I randomly initialized a population of 50 particles, some of which are already close to the solution.

Now let’s implement the PSO algorithm itself. I commented each row in the code, but if you have any questions, feel free to ask in the comments below.

"""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}')

After 30 iteration we’ve got this:

As you can see the algorithm fell into the local minimum, which is not what we wanted. That’s why we need to tune our hypoparameters and start again. This time I decided to set inertia weight **w=0.8**, thus, now the previous velocity has a greater impact on the current state.

And voila, we reached the global minimum of the function. I strongly encourage you to play around with POP_SIZE, c₁ and c₂. It’ll allow you to gain a better understanding of the code and the idea behind PSO. If you’re interested you can complicate the task and optimize some 3D function and make a nice visualization.

===========================================

[1]Shi Y. Particle swarm optimization //IEEE connections. — 2004. — Т. 2. — №. 1. — С. 8–13.

===========================================

*All my articles on Medium are free and open-access, that’s why I’d really appreciate if you followed me here!*

P.s. I’m extremely passionate about (Geo)Data Science, ML/AI and Climate Change. So if you want to work together on some project pls contact me in LinkedIn.

🛰️Follow for more🛰️

Particle Swarm Optimization (PSO) from scratch. Simplest explanation in python was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.

**Denial of responsibility!**Techno Blender is an automatic aggregator of the all world’s media. In each content, the hyperlink to the primary source is specified. All trademarks belong to their rightful owners, all materials to their authors. If you are the owner of the content and do not want us to publish your materials, please contact us by email – [email protected]. The content will be deleted within 24 hours.