Particle Swarm Optimization Basics

The implementation presented here is the original PSO algorithm as presented in [Poli2007]. From Wikipedia definition of PSO

PSO optimizes a problem by having a population of candidate solutions, here dubbed particles, and moving these particles around in the search-space according to simple mathematical formulae. The movements of the particles are guided by the best found positions in the search-space which are updated as better positions are found by the particles.

Modules

Before writing functions and algorithms, we need to import some module from the standard library and from DEAP.

import operator
import random
from itertools import imap

from deap import base
from deap import benchmarks
from deap import creator
from deap import tools

Representation

The particle’s goal is to maximize the return value of the fonction at its position.

PSO particles are essentially described as a positions in a search-space of D dimensions. Each particle also has a vector representing the speed of the particle in each dimension. Finally, each particle keeps a reference to the best state in which it has been so far.

This translates in DEAP by the following two lines of code :

creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Particle", list, fitness=creator.FitnessMax, speed=list, 
    smin=None, smax=None, best=None)

Here we create two new objects in the creator space. First, we create a FitnessMax object, and we specify the weights to be (1.0,), this means we want to maximise the value of the fitness of our particles. The second object we create represent our particle. We defined it as a list to which we add five attributes. The first attribute is the fitness of the particle, the second is the speed of the particle which is also goind to be a list, the third and fourth are the limit of the speed value, and the fifth attribute will be a reference to a copy of the best state the particle has been so far. Since the particle has no final state until at has been evaluated, the reference is set to None. The speed limits are also set to None to allow configuration via the function generate() presented in the next section.

Operators

PSO original algorithm use three operators : initializer, updater and evaluator. The initialization consist in generating a random position and a random speed for a particle. The next function create a particle and initialize its attributes, except for the attribute best, which will be set only after evaluation :

def generate(size, pmin, pmax, smin, smax):
    part = creator.Particle(random.uniform(pmin, pmax) for _ in xrange(size)) 
    part.speed = [random.uniform(smin, smax) for _ in xrange(size)]
    part.smin = smin
    part.smax = smax
    return part

The function updateParticle() first computes the speed, then limits the speed values between smin and smax, and finally computes the new particle position.

def updateParticle(part, best, phi1, phi2):
    u1 = (random.uniform(0, phi1) for _ in range(len(part)))
    u2 = (random.uniform(0, phi2) for _ in range(len(part)))
    v_u1 = imap(operator.mul, u1, imap(operator.sub, part.best, part))
    v_u2 = imap(operator.mul, u2, imap(operator.sub, best, part))
    part.speed = map(operator.add, part.speed, imap(operator.add, v_u1, v_u2))
    for i, speed in enumerate(part.speed):
        if speed < part.smin:
            part.speed[i] = part.smin
        elif speed > part.smax:
            part.speed[i] = part.smax
    part[:] = map(operator.add, part, part.speed)

The operators are registered in the toolbox with their parameters. The particle value at the beginning are in the range [-100, 100] (pmin and pmax), and the speed is limited in the range [-50, 50] through all the evolution.

The evaluation function h1() is from [Knoek2003]. The function is already defined in the benchmarks module, so we can register it directly.

toolbox = base.Toolbox()
toolbox.register("particle", generate, size=2, pmin=-6, pmax=6, smin=-3, smax=3)
toolbox.register("population", tools.initRepeat, list, toolbox.particle)
toolbox.register("update", updateParticle, phi1=2.0, phi2=2.0)
toolbox.register("evaluate", benchmarks.h1)

Algorithm

Once the operators are registered in the toolbox, we can fire up the algorithm by firstly creating a new population, and then apply the original PSO algorithm. The variable best contains the best particle ever found (it known as gbest in the original algorithm).

def main():
    pop = toolbox.population(n=5)
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("Avg", tools.mean)
    stats.register("Std", tools.std)
    stats.register("Min", min)
    stats.register("Max", max)
    
    logger = tools.EvolutionLogger(["gen", "evals"] + stats.functions.keys())
    logger.logHeader()

    GEN = 1000
    best = None

    for g in xrange(GEN):
        for part in pop:
            part.fitness.values = toolbox.evaluate(part)
            if not part.best or part.best.fitness < part.fitness:
                part.best = creator.Particle(part)
                part.best.fitness.values = part.fitness.values
            if not best or best.fitness < part.fitness:
                best = creator.Particle(part)
                best.fitness.values = part.fitness.values
        for part in pop:
            toolbox.update(part, best)

        # Gather all the fitnesses in one list and print the stats
        stats.update(pop)
        logger.logGeneration(gen=g, evals=len(pop), stats=stats)
    
    return pop, stats, best

Conclusion

The full PSO basic example can be found

This is a video of the algorithm in action, plotted with matplotlib. The red dot represents the best solution found so far.

References

[Poli2007]

Ricardo Poli, James Kennedy and Tim Blackwell, “Particle swarm optimization an overview”. Swarm Intelligence. 2007; 1: 33–57

[Knoek2003]

Arthur J. Knoek van Soest and L. J. R. Richard Casius, “The merits of a parallel genetic algorithm in solving hard optimization problems”. Journal of Biomechanical Engineering. 2003; 125: 141–146