Particle Swarm Minimiser

This is an implementation of a particle swarm optimisation algorithm, which will find the minimum of two functions as well as report the corresponding input values. The two functions to minimise were

$f(x, y) = 1 \,+ \,&(-13 + x - y^3 + 4y^2 - 2y)^2 \,+ \ &(-29 + x + y^3 + y^2 - 14y)^2$

and

$f=- &\begin{bmatrix}15 & 27 & 36 & 18 & 12\end{bmatrix} \mathbf{x} \, +\ \mathbf{x} ^ \top{} &\begin{bmatrix*}[r] 35 & -20 & -10 & 32 & -10\ -20 & 40 & -6 & -31 & 32\ -10 & -6 & 11 & -6 & -10\ 32 & -31 & -6 & 38 & -20\ -10 & 32 & -10 & -20 & 31 \end{bmatrix*} \mathbf{x}\,$

The algorithm is inspired by how for example bird flocks or fish schools behave. They are an interesting study, as there are no apparent leaders that the other members follow, yet the swarm behaves coherent.

In particle swarm optimisation, a common model is the Boids model (the name "Boids" is taken from "bird-like objects"), which is the model that was implemented in this program. In it, we let $\mathbf{S}\,$ denote a swarm of boids, hereafter called particles:

$\mathbf{S} = \{p_i | i = 1, ..., N \}\,$ .

In the model, each particle p_i only has a limited visual range, and is unable to perceive swarm mates that are far away. The so called visibility sphere $\mathbf{V}_i$ of particle is introduced as

$\mathbf{V}_i = \{ p_j : || \mathbf{x}_j - \mathbf{x}_i || < r, \:\: j \neq i \}\,$ ,

where is a constant. Each particle p_i has a position $\mathbf{x}_i$ , a velocity $\mathbf{v}_i$ and an acceleration $\mathbf{a}_i$ , and in each iteration of the algorithm, they are updated using standard Euler integration:

$\mathbf{v}_i &\leftarrow \mathbf{v}_i + \mathbf{a}_i \Delta t\ \mathbf{x}_i &\leftarrow \mathbf{x}_i + \mathbf{v}_i \Delta t$ ,

where $\Delta t$ is the time step (here set to 1). Each particle is influenced by three movement tendencies that together determine its acceleration. These tendencies are cohesion, alignment and separation.

Cohesion
This is the tendency of a particle to stay near the centre of the swarm. Let $\rho_i$ denote the centre of density of the visibility sphere of particle . If there are particles in $\mathbf{V}_i$ , then
$\rho_i = \frac{1}{k_i} \sum_{p_j \in \mathbf{V}_i} \mathbf{x}_i$
The steering vector representing cohesion is defined as
$\mathbf{c}_i = \frac{1}{T^2}\left( \rho_i - \mathbf{x}_i \right)$
where is a time constant. If there are no particles in $\mathbf{V}_i$ , then $\mathbf{c}_i = \mathbf{0}\,$
Alignment
This is the tendency of particles to align their velocities with those of their nearby swarm mates. It is defined as
$\mathbf{l}_i = \frac{1}{Tk_i} \sum_{p_j \in \mathbf{V}_i} \mathbf{v}_j$
where is a time constant and are the number of particles in $\mathbf{V}_i$ . If there are no particles in $\mathbf{V}_i$ , then $\mathbf{l}_i = \mathbf{0}\,$ .
Separation
Separation is needed in order to avoid collision with nearby swarm mates, and is defined as
$\mathbf{s}_i = \frac{1}{T^2} \sum_{p_j \in \mathbf{V}_i} (\mathbf{x}_i - \mathbf{x}_j)$
where once again is a time constant. If there are no particles in $\mathbf{V}_i$ , then $\mathbf{s}_i = \mathbf{0}\,$ .

The acceleration of particle is obtained by combining these steering vectors, by the following equation:
$\mathbf{a}_i = C_c\mathbf{c}_i + C_l\mathbf{l}_i + C_s\mathbf{s}_i$ ,
where C_c , C_l and C_s are constants between 0 and 1 that control the relative impact of the steering vectors.

The algorithm works as follows:

Initialisation
Initialise the positions and velocities of each particle. Let denote the swarm size and the number of variables in the problem to solve. The positions are randomly initialised, using uniform sampling in the range between $x_{\textrm{min} }\,$ and $x_{\textrm{max} }\,$ .

Initialise positions of all the particles as follows:
$x_{ij} = x_{\textrm{min} } + r(x_{\textrm{max} } - x_{\textrm{min} })$
where $x_{ij}\,$ denotes component of the position of particle and is a random number in the range 0 to 1.

Initialise the velocities of all the particles as follows:
$v_{ij} = \frac{\alpha}{\Delta t} \left ( - \frac{x_{\textrm{max} } - x_{\textrm{min} } }{2} + r \left ( x_{\textrm{max} } - x_{\textrm{min} } \right ) \right )$
where $v_{ij}\,$ denotes component of the velocity of particle , is a random number between 0 and 1, $\alpha$ is a constant between 0 and 1, and $\Delta t$ is the time step length. In the implementation, $\alpha$ and $\Delta t$ are both set to 1.
Evaluation
The evaluation part consists of computing the objective function value for each $\mathbf{x}_i$ . In other words, input $\mathbf{x}_i$ into the objective function, in this case either
$f(x, y) = 1 \,+ \,&(-13 + x - y^3 + 4y^2 - 2y)^2 \,+ \ &(-29 + x + y^3 + y^2 - 14y)^2$
or
$f=- &\begin{bmatrix}15 & 27 & 36 & 18 & 12\end{bmatrix} \mathbf{x} \, +\ \mathbf{x} ^ \top{} &\begin{bmatrix*}[r] 35 & -20 & -10 & 32 & -10\ -20 & 40 & -6 & -31 & 32\ -10 & -6 & 11 & -6 & -10\ 32 & -31 & -6 & 38 & -20\ -10 & 32 & -10 & -20 & 31 \end{bmatrix*} \mathbf{x}\,$
Update the best positions
The objective of the algorithm is to reach optimal values of the objective function. Thus, the algorithm must keep track of the particles' performance. Both the best position $\mathbf{x}_{i}^{\textrm{pb} }\,$ of particle , and the best performance $\textbf{x}_{\textrm{sb} }\,$ of any particle in the swarm is stored. The former means comparing the position of particle in the current iteration to how well it has previously performed, and updating the record if needed. The latter can be carried out in different ways, but in this implementation, $\textbf{x}_{\textrm{sb} }\,$ is the best-ever position of the whole swarm. The following equations capture this:
$\textrm{if}\:\: & f(\mathbf{x}_i) < f(\mathbf{x}_{i}^{\textrm{pb} }) \:\: \textrm{then} \:\: \mathbf{x}_{i}^{\textrm{pb} } \leftarrow \mathbf{x}_i\ \textrm{if}\:\: & f(\mathbf{x}_i) < f(\mathbf{x}^ {\textrm{sb} })\,\:\:\textrm{then} \:\: \mathbf{x} ^{\textrm{sb} }\,\leftarrow \mathbf{x}_i$
Update particle velocities and positions
Update the velocity of particle as follows, where goes from 1 to :
$v_{ij} \leftarrow wv_{ij} + c_1 q \left ( \frac {x_{ij}^\textrm{pb} - x_{ij} }{\Delta t} \right ) + c_2 r \left ( \frac {x_{j}^\textrm{sb} - x_{ij} }{\Delta t} \right )$
where and are uniform random numbers between 0 and 1, is the so called inertia weight and and are positive constants. The -term is sometimes called the cognitive component, and measures the degree of self-confidence of a particle, i.e. how much it trusts its own previous performance in order to get a better result. The -term is sometimes called the social component, and measures how much a particle trusts others to find better candidate solutions.

The inertia weight handles the trade off between exploration () and exploitation (). Initially, , and after each iteration, is modified as follows:
$w \leftarrow \beta w$ ,
where $\beta$ is 0.99. For every iteration, is decreased in this manner, until it reaches a minimum value of 0.35, where it will be kept fixed. This procedure makes the algorithm explorative in the early stages, and later on turns more and more into exploitation.

If the velocities are allowed to grow with no boundary, the swarm will cease to be coherent and will expand indefinitely. Therefore, a maximal velocity $v_{\textrm{max} }\,$ is introduced, restricting velocities by
$v_{ij} = \begin{cases} v_{ij} & \text{if $|v_{ij}| < v_{\textrm{max} }$}\ \textrm{sign}(v_{ij})v_{\textrm{max} } & \text{otherwise} \end{cases}\,$
where the sign-function is used to maintain the direction, and is defined as
$\textrm{sign}(a) = \begin{cases} -1 & \textrm{if}\:\: a < 0\ 0 & \textrm{if}\:\: a = 0\ 1 & \textrm{if}\:\: a > 0 \end{cases}\,$
The position of particle is then updated as follows, where goes from 1 to :
$x_{ij} \leftarrow x_{ij} + v_{ij} \Delta t$
Repeat
Repeat from step 2, until the termination criterion has been reached.

The code is available at my GitHub account.

Screenshot of Matlab running the program:

Wikimedia contributions

Particle Swarm Minimiser