Particle Swarm Optimisation

Published

February 13, 2024

Particle swarm optimisation is inspired by the social behaviour of birds and fish. It is a population-based optimisation algorithm that is used to find the best solution to a problem. The algorithm is based on the idea of a swarm of particles that move around in a search space, and each particle represents a potential solution to the problem. The particles move around the search space, and their movement is influenced by their own experience and the experience of the other particles in the swarm.

The algorithm will be explained for the following unconstrained optimisation problem:

\[ \min_{\mathbf{x}\in \mathbb{R}^n} f(\mathbf{x}) \]

The pseudocode for the algorithm is as follows:

  1. Initialise the swarm of particles with random positions and velocities \(\mathbf{x}_i\) and \(\mathbf{v}_i\) where \(i \in [1, n_p]\) and \(n_p\) is the number of particles in the swarm.

  2. Evaluate the objective function for each particle (\(f(\mathbf{x}_i)\)) and update the best position (\(\mathbf{p}_i\)) and best objective function value (\(f(\mathbf{p}_i)\)) for each particle.

  3. Update the best position and best objective function value for the swarm.

  4. Update the velocity and position of each particle using the following equations:

\[ \mathbf{v}_{i+1} = w\,\mathbf{v}_i + c_1\,r_1\,(\mathbf{p}_i - \mathbf{x}_i) + c_2\,r_2(\mathbf{g}- \mathbf{x}_i) \mathbf{x}_{i+1} = \mathbf{x}_i + \mathbf{v}_{i+1} \]

where \(w\) is the inertia weight, \(c_1\) and \(c_2\) are the acceleration coefficients, and \(r_1\) and \(r_2\) are random numbers between 0 and 1.

  1. Repeat steps 2-4 until a stopping criterion is met.

The algorithm is sensitive to the choice of parameters, and the performance of the algorithm can be improved by tuning the parameters. The inertia weight \(w\) controls the influence of the previous velocity on the current velocity, and the acceleration coefficients \(c_1\) and \(c_2\) control the influence of the best position of the particle and the best position of the swarm on the velocity of the particle.

Best position and best objective function value for each particle are updated as follows:

\[ \mathbf{p}_i = \begin{cases} \mathbf{x}_i & \text{if } f(\mathbf{x}_i) < f(\mathbf{p}_i) \\ \mathbf{p}_i & \text{otherwise} \end{cases} \]

\[ f(\mathbf{p}_i) = \min(f(\mathbf{x}_i), f(\mathbf{p}_i)) \]

The best position and best objective function value for the swarm are updated as follows:

\[ \mathbf{g}= \min_{i \in [1,n_p]} f(\mathbf{p}_i) \]

\[ f(\mathbf{g}) = \min_{i \in [1,n_p]} f(\mathbf{p}_i) \]

The stopping criterion can be a maximum number of iterations, a maximum number of function evaluations, or a minimum change in the objective function value.

Here is an example python code that demonstrates the particle swarm optimisation algorithm for a bi-model objective function. The initial popolation of \(10\) is generated randomly and the algorithm is run for \(100\) iterations. The inertia weight is set to \(0.5\), and the acceleration coefficients are set to \(1.5\). The objective function is defined as:

\[ f(\mathbf{x}) = \sin(x_1) + \cos(x_2) \]

The objective function has two global minima at \(\mathbf{x}= (n\pi, m\pi)\) where \(n, m \in \mathbb{Z}\). The algorithm is run for \(100\) iterations, and the best objective function value is plotted as a function of the iteration number.

The output of the code is an image that shows the evolution of the population through the iterations. For the sake of clarity, the populations are plotted after every \(10\) iterations. The best objective function value is plotted as a function of the iteration number.

Code 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

# Define the Rosenbrock function
def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# Define the PSO algorithm
def particle_swarm_optimisation(f, swarm, num_iterations, w=0.5, c1=0.05, c2=0.05):
    positions = []
    velocities = np.random.uniform(-1, 1, swarm.shape)
    personal_best_positions = np.copy(swarm)
    personal_best_values = np.array([f(x) for x in swarm])
    global_best_position = swarm[np.argmin(personal_best_values)]

    for i in range(num_iterations):
        for j in range(len(swarm)):
            # Update velocity
            velocities[j] = (w * velocities[j] +
                             c1 * np.random.random() * (personal_best_positions[j] - swarm[j]) +
                             c2 * np.random.random() * (global_best_position - swarm[j]))

            # Update position
            swarm[j] += velocities[j]

            # Update personal best
            if f(swarm[j]) < personal_best_values[j]:
                personal_best_positions[j] = swarm[j]
                personal_best_values[j] = f(swarm[j])

            # Update global best
            if f(swarm[j]) < f(global_best_position):
                global_best_position = swarm[j]

        if i % 10 == 0:
            positions.append(np.copy(swarm))

    return global_best_position, np.array(positions)

# Initialize the swarm
np.random.seed(0)
swarm = np.random.uniform(-3, 3, (50, 2))

# Run the PSO algorithm
best_solution, positions = particle_swarm_optimisation(rosenbrock, swarm, 100)

# Create a scatter plot for the initial particle positions
fig, ax = plt.subplots()
scat = ax.scatter(positions[0][:, 0], positions[0][:, 1])

# Update function for the animation
def update(i):
    # Update the scatter plot for the current particle positions
    scat.set_offsets(positions[i])

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=len(positions), interval=1)
plt.close()
HTML(ani.to_jshtml())

# # Display the animation
# plt.show()
ModuleNotFoundError: No module named 'numpy'
Code 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Define the Rosenbrock function
def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# Define the PSO algorithm
def particle_swarm_optimisation(f, swarm, num_iterations, w=0.5, c1=1, c2=1):
    positions = []
    velocities = np.random.uniform(-1, 1, swarm.shape)
    personal_best_positions = np.copy(swarm)
    personal_best_values = np.array([f(x) for x in swarm])
    global_best_position = swarm[np.argmin(personal_best_values)]

    for i in range(num_iterations):
        for j in range(len(swarm)):
            # Update velocity
            velocities[j] = (w * velocities[j] +
                             c1 * np.random.random() * (personal_best_positions[j] - swarm[j]) +
                             c2 * np.random.random() * (global_best_position - swarm[j]))

            # Update position
            swarm[j] += velocities[j]

            # Update personal best
            if f(swarm[j]) < personal_best_values[j]:
                personal_best_positions[j] = swarm[j]
                personal_best_values[j] = f(swarm[j])

            # Update global best
            if f(swarm[j]) < f(global_best_position):
                global_best_position = swarm[j]

        if i % 10 == 0:
            positions.append(np.copy(swarm))

    return global_best_position, np.array(positions)

# Initialize the swarm
np.random.seed(0)
swarm = np.random.uniform(-3, 3, (50, 2))

# Run the PSO algorithm
best_solution, positions = particle_swarm_optimisation(rosenbrock, swarm, 1000)

# Create a scatter plot for the initial particle positions
fig, ax = plt.subplots()

# Create a grid of points
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)

# Compute the function value at each point
Z = rosenbrock([X, Y])

# Create a contour plot
contour = ax.contour(X, Y, Z, levels=np.logspace(0, 5, 35), cmap='jet')
scat = ax.scatter(positions[0][:, 0], positions[0][:, 1], color='k')

# Update function for the animation
def update(i):
    # Update the scatter plot for the current particle positions
    scat.set_offsets(positions[i])

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=len(positions), interval=200)

# Save the animation as a GIF file
ani.save('animation.gif', writer='pillow')

# Display the animation
plt.show()
ModuleNotFoundError: No module named 'numpy'

Animation