PSO: Particle Swarm Optimization

Particle Swarm Optimization was proposed in 1995 by Kennedy and Eberhart [18] based on the simulating of social behavior. The algorithm uses a swarm of particles to guide its search. Each particle has a velocity and is influenced by locally and globally best-found solutions. Many different implementations have been proposed in the past and, therefore, it is quite difficult to refer to THE correct implementation of PSO. However, the general concepts shall be explained in the following.

Given the following variables:

  • \(X_{d}^{(i)}\) d-th coordinate of i-th particle’s position

  • \(V_{d}^{(i)}\) d-th coordinate of i-th particle’s velocity

  • \(\omega\) Inertia weight

  • \(P_{d}^{(i)}\) d-th coordinate of i-th particle’s personal best

  • \(G_{d}^{(i)}\) d-th coordinate of the globally (sometimes also only locally) best solution found

  • \(c_1\) and \(c_2\) Two weight values to balance exploiting the particle’s best \(P_{d}^{(i)}\) and swarm’s best \(G_{d}^{(i)}\)

  • \(r_1\) and \(r_2\) Two random values being create for the velocity update

The velocity update is given by:

\begin{equation} V_{d}^{(i)} = \omega \, V_{d}^{(i)} \;+\; c_1 \, r_1 \, \left(P_{d}^{(i)} - X_{d}^{(i)}\right) \;+\; c_2 \, r_2 \, \left(G_{d}^{(i)} - X_{d}^{(i)}\right) \end{equation}

The corresponding position value is then updated by:

\begin{equation} X_{d}^{(i)} = X_{d}^{(i)} + V_{d}^{(i)} \end{equation}

The social behavior is incorporated by using the globally (or locally) best-found solution in the swarm for the velocity update. Besides the social behavior, the swarm’s cognitive behavior is determined by the particle’s personal best solution found. The cognitive and social components need to be well balanced to ensure that the algorithm performs well on a variety of optimization problems. Thus, some effort has been made to determine suitable values for \(c_1\) and \(c_2\). In pymoo both values are updated as proposed in [19]. Our implementation deviates in some implementation details (e.g. fuzzy state change) but follows the general principles proposed in the paper.

[1]:
from pymoo.algorithms.soo.nonconvex.pso import PSO, PSOAnimation
from pymoo.factory import Ackley
from pymoo.optimize import minimize

problem = Ackley()

algorithm = PSO(max_velocity_rate=0.025)

res = minimize(problem,
               algorithm,
               callback=PSOAnimation(fname="pso.mp4", nth_gen=5),
               seed=1,
               save_history=True,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))
Best solution found:
X = [-1.49794464e-06  3.25336775e-06]
F = [1.01307867e-05]

Assuming you have our third-party library pyrecorder installed you can create an animation for a two-dimensional problem easily. Below we provide the code to observe the behavior of the swarm. We have reduced to set the maximum velocity to max_velocity_rate=0.025 for illustration purposes. Otherwise, the algorithm converges even quicker to the global optimum of the Ackley function with two variables.

In general, the PSO algorithm can be used by execute the following code. For the available parameters please see the API description below.

[2]:
from pymoo.algorithms.soo.nonconvex.pso import PSO
from pymoo.factory import Rastrigin
from pymoo.optimize import minimize

problem = Rastrigin()

algorithm = PSO()

res = minimize(problem,
               algorithm,
               seed=1,
               verbose=False)

print("Best solution found: \nX = %s\nF = %s" % (res.X, res.F))
Best solution found:
X = [-8.84624648e-07  1.71712265e-06]
F = [7.40215e-10]

API

class pymoo.algorithms.soo.nonconvex.pso.PSO(self, pop_size=25, sampling=LHS(), w=0.9, c1=2.0, c2=2.0, adaptive=True, initial_velocity='random', max_velocity_rate=0.20, pertube_best=True, repair=NoRepair(), display=PSODisplay(), **kwargs)
Parameters
pop_sizeThe size of the swarm being used.
samplingSampling, Population, numpy.array

The sampling process defines the initial set of solutions which are the starting point of the optimization algorithm. Here, you have three different options by passing

(i) A Sampling implementation which is an implementation of a random sampling method.

(ii) A Population object containing the variables to be evaluated initially OR already evaluated solutions (F needs to be set in this case).

(iii) Pass a two dimensional numpy.array with (n_individuals, n_var) which contains the variable space values for each individual.

adaptivebool

Whether w, c1, and c2 are changed dynamically over time. The update uses the spread from the global optimum to determine suitable values.

wfloat

The inertia F to be used in each iteration for the velocity update. This can be interpreted as the momentum term regarding the velocity. If adaptive=True this is only the initially used value.

c1float

The cognitive impact (personal best) during the velocity update. If adaptive=True this is only the initially used value.

c2float

The social impact (global best) during the velocity update. If adaptive=True this is only the initially used value.

initial_velocitystr - (‘random’, or ‘zero’)

How the initial velocity of each particle should be assigned. Either ‘random’ which creates a random velocity vector or ‘zero’ which makes the particles start to find the direction through the velocity update equation.

max_velocity_ratefloat

The maximum velocity rate. It is determined variable (and not vector) wise. We consider the rate here since the value is normalized regarding the xl and xu defined in the problem.

pertube_bestbool

Some studies have proposed to mutate the global best because it has been found to converge better. Which means the population size is reduced by one particle and one function evaluation is spend additionally to permute the best found solution so far.