Particle Swarm Optimization in MATLAB

Particle Swarm Optimization (PSO) searches for a low-cost solution with a population of moving candidates called particles. Each particle combines its current velocity with directions toward its own best position and the best position found by the swarm.

Kennedy and Eberhart introduced PSO in 1995, drawing on simplified models of collective behavior. The MATLAB implementation here turns that idea into a complete optimization loop, with a Sphere test problem, position bounds, velocity limits, and a convergence plot. For background, see the Wikipedia article on PSO (opens in a new tab).

How PSO works

Each particle $i$ has a position $x_i$, which is a candidate solution, and a velocity $v_i$. It remembers the best position it has found so far, its personal best $p_i$, and the swarm shares the best position found by any particle, the global best $g$. In every iteration $t$, each particle first updates its velocity and then moves:

Each term of the velocity update has a clear role:

• $w\,v_i^{t}$ is the inertia term. It keeps the particle moving in its current direction, and the inertia weight $w$ sets how strongly.
• $c_1 r_1 \odot (p_i - x_i^{t})$ is the cognitive term. It pulls the particle back toward its own best position.
• $c_2 r_2 \odot (g - x_i^{t})$ is the social term. It pulls the particle toward the best position of the swarm.

Here $c_1$ and $c_2$ are the personal and social acceleration coefficients, and $r_1$ and $r_2$ are vectors of random numbers drawn uniformly from $[0, 1]$ in every iteration. The symbol $\odot$ means element by element multiplication, so each dimension gets its own random weight.

After the move, each particle evaluates the cost function $f$ and updates its memory. For a minimization problem, the personal best changes only when the new position is better:

The global best $g$ is then the best of all personal bests.

Inertia weight

A large inertia weight favors exploration of the search space, and a small one favors exploitation near the best solutions found so far. A common strategy is to start with $w$ close to 1 and reduce it a little in every iteration with a damping ratio $w_{\text{damp}}$ slightly below 1, for example 0.99:

Constriction coefficients

The constriction formulation of Clerc and Kennedy (2002) gives a standard way to choose the velocity coefficients. Pick two positive numbers $\varphi_1$ and $\varphi_2$ with $\varphi = \varphi_1 + \varphi_2 > 4$, and compute the constriction factor

The coefficients are then $w = \chi$, $c_1 = \chi\varphi_1$ and $c_2 = \chi\varphi_2$. The usual choice $\varphi_1 = \varphi_2 = 2.05$ gives $\varphi = 4.1$, $\chi \approx 0.7298$ and $c_1 = c_2 \approx 1.4962$.

Implementation in MATLAB

This post shares a complete implementation of PSO in MATLAB. The code is structured, so if you know MATLAB, you should find it easy to use in your research projects.

What is in the download

The PSO folder holds three files:

• pso.m: the PSO script, with the problem definition, the parameters, the initialization, the main loop and the results.
• Sphere.m: the sample cost function, $f(x) = \sum_{i=1}^{n} x_i^2$, whose minimum is at the origin.
• main.m: a one line entry point that runs pso.m.

To run the code, open the PSO folder in MATLAB and run main.m. The program prints the best cost of every iteration and then draws the convergence plot (best cost against iteration) on a logarithmic scale with semilogy.

The script creates the first positions with unifrnd, a function of the Statistics and Machine Learning Toolbox, so it probably needs that toolbox.

Parameters in the code

The parameters are set at the top of pso.m:

To use the constriction coefficients instead, uncomment the block in pso.m that sets $\varphi_1 = \varphi_2 = 2.05$, computes $\chi$, and sets $w = \chi$, $c_1 = \chi\varphi_1$, $c_2 = \chi\varphi_2$ and wdamp to 1, and comment out the four lines above it.

Bounds and velocity limits

Two details of the code go beyond the basic equations:

• The velocity of each particle is clipped to $[-v_{\max}, v_{\max}]$, with $v_{\max} = 0.1\,(\text{VarMax} - \text{VarMin})$.
• After a move, any position outside the bounds is set back to the nearest bound, and the velocity in that dimension is reversed, so the particle moves back into the search space.

To solve your own problem, write a cost function that takes a row vector and returns a number, and set CostFunction, nVar, VarMin and VarMax in the problem definition section of pso.m.

Check the result

Start with the supplied Sphere problem. Its minimum cost is 0 at the origin, which lies inside the default bounds. Check the final position and cost together, then repeat the run to see the effect of random initialization. When changing the problem, verify that the cost function accepts a row vector, returns one finite scalar, and evaluates the intended units and constraints. Constriction coefficients and a small final cost do not establish global optimality for an arbitrary objective.

Applications

PSO needs only the values of the cost function, not its gradient, so it works on problems where the cost is noisy or has many local minima. Typical uses are the tuning of controllers and model parameters, engineering design, and the training of neural networks and other models. Yarpiz has several application posts that use PSO, for example optimal inventory control, economic dispatching and robot path planning.

Related posts

• Particle Swarm Optimization (PSO) in Python is a Python version of the same code.
• YPEA is a MATLAB toolbox with PSO and 14 other algorithms that share one problem definition.
• MOPSO in MATLAB is a MATLAB implementation of Multi-Objective Particle Swarm Optimization.

Mostapha Kalami Heris, Particle Swarm Optimization in MATLAB (YPEA102), Yarpiz, 2015. https://yarpiz.com/particle-swarm-optimization-ypea102

@misc{yarpiz-ypea102,
  author = {{Mostapha Kalami Heris}},
  title = {{Particle Swarm Optimization in MATLAB}},
  howpublished = {Yarpiz},
  year = {2015},
  note = {Code YPEA102},
  url = {https://yarpiz.com/particle-swarm-optimization-ypea102}
}