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.
Before writing functions and algorithms, we need to import some module from the standard library and from DEAP.
import operator
import random
from deap import base
from deap import benchmarks
from deap import creator
from deap import tools
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.
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 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 speed < part.smin:
part.speed[i] = part.smin
elif speed > part.smax:
part.speed[i] = part.smax
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)
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)
column_names = ["gen", "evals"]
column_names.extend(stats.functions.keys())
logger = tools.EvolutionLogger(column_names)
logger.logHeader()
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
stats.update(pop)
logger.logGeneration(gen=g, evals=len(pop), stats=stats)
return pop, stats, best
The full PSO basic example can be found here : [source code].
This is a video of the algorithm in action, plotted with matplotlib. The red dot represents the best solution found so far.
[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 |