# 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

    import numpy
    import math

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

## Representation

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

PSO particles are essentially described as 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 going 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 it 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 uses three operators : initializer, updater and evaluator. The initialization consist in generating a random position and a random speed for a particle. The next function creates a particle and initializes 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 range(size))
        part.speed = [random.uniform(smin, smax) for _ in range(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 = map(operator.mul, u1, map(operator.sub, part.best, part))
        v_u2 = map(operator.mul, u2, map(operator.sub, best, part))
        part.speed = list(map(operator.add, part.speed, map(operator.add, v_u1, v_u2)))
        for i, speed in enumerate(part.speed):
            if abs(speed) < part.smin:
                part.speed[i] = math.copysign(part.smin, speed)
            elif abs(speed) > part.smax:
                part.speed[i] = math.copysign(part.smax, speed)
        part[:] = list(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 is 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", numpy.mean)
        stats.register("std", numpy.std)
        stats.register("min", numpy.min)
        stats.register("max", numpy.max)

        logbook = tools.Logbook()
        logbook.header = ["gen", "evals"] + stats.fields

        GEN = 1000
        best = None

        for g in range(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
            logbook.record(gen=g, evals=len(pop), **stats.compile(pop))
            print(logbook.stream)

        return pop, logbook, best

## Conclusion

The full PSO basic example can be found here : examples/pso/basic.py.

## 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