python DEAP PSO 算法的学习

翻译自：http://deap.readthedocs.io/en/master/examples/pso_basic.html

粒子群优化基础

这里的实现是基本PSO算法的实现。wikipedia的PSO的定义如下：
PSO通过拥有一个种群规模的候选解决方案来优化一个问题，这里的候选解决方案指的是粒子，在解空间中根据简单的数学公式移动这些粒子。粒子的移动被解空间的最优位置引导着，如果粒子发现更好的位置就会更新最优位置。

模块

在写函数和算法之前，需要从标准库和DEAP中导入一些模块。

import operator
import random

import numpy

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

实现

粒子的目标是在当前位置最大化函数的返回值。
PSO算法中的粒子实质上被描述为D维搜索空间上的一个位置。每一个粒子用一个向量来表示粒子在每一维的速度。最终，每个粒子参照目前为止最好的状态进行移动。
在DEAP中转化为如下两行代码：

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

在creator创建了两个新的对象。首先，创建一个FitnessMax对象，指定weights为（1.0，）,这意味着我们想要最大化粒子的fitness的值（因为是正数）。第二个对象表示粒子，我们定义这个例子为一个list，并添加5个属性。第一个属性是粒子的fitness适应度，第二个是粒子的速度，也是一个list，第三个和第四个是速度的临界值，第五个是粒子迄今为止最佳状态的副本的引用。因为粒子直到被评估之前都没有最终的状态，所以引用值初始设置为None。速度的临界值也被设置为None，可以通过函数generate()设置这个值，下一部分将会介绍。

operators 算子

PSO基本算法只用三个算子：初始化，更新，评估。初始化是为每一个粒子产生一个随机的位置和随机的速度。接下来这个函数创建一个粒子并初始化属性，除了属性best，这个需要在评估以后才能设置。

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

函数updateParticle()先计算速度，然后限制速度在smin和smax之间，最后计算新的粒子位置。

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 par

这些算子使用参数在toolbox注册。粒子值定义在范围 [-100, 100] (pmin 和pmax)，速度限制在范围 [-50, 50] 。评估函数 h1() 来自 [Knoek2003]论文。这个函数已经在benchmarks模块中定义了，所以可以直接注册并使用。

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)

算法

一旦算子在toolbox中注册以后，我们就可以使用算法来创建一个新的种群，然后应用原始PSO算法。变量best包含迄今为止发现的最优的粒子（在原始算法中表示gbest全局最优解）。

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

结论

完整的PSO例子在这里：https://github.com/DEAP/deap/blob/09b2562aad70c6f171ebca2bdac0c30387f5e8f5/examples/pso/basic.py
链接的代码如下：

#    This file is part of DEAP.
#
#    DEAP is free software: you can redistribute it and/or modify
#    it under the terms of the GNU Lesser General Public License as
#    published by the Free Software Foundation, either version 3 of
#    the License, or (at your option) any later version.
#
#    DEAP is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#    GNU Lesser General Public License for more details.
#
#    You should have received a copy of the GNU Lesser General Public
#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.

import operator
import random

import numpy

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

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

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

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

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)

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

if __name__ == "__main__":
    main()