粒子群优化基础

这里介绍的实现是原始的PSO算法,如[Poli2007]_中所述。根据维基百科对PSO的定义

PSO 通过拥有一组候选解(这里称为粒子)来优化问题,并根据简单的数学公式在搜索空间中移动这些粒子。粒子的移动由搜索空间中找到的最佳位置引导,这些最佳位置会随着粒子找到更好的位置而更新。

模块

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

import operator
import random

import numpy
import math

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

表示

粒子的目标是最大化其位置上的函数返回值。

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,),这意味着我们希望最大化我们粒子的适应度值。我们创建的第二个对象代表我们的粒子。我们将其定义为一个 list,并向其中添加了五个属性。第一个属性是粒子的适应度,第二个是粒子的速度,它也将是一个列表,第三和第四是速度值的限制,第五个属性将是对粒子迄今为止所处的最佳状态的副本的引用。由于粒子在被评估之前没有最终状态,因此引用被设置为 None。速度限制也设置为 None,以便通过下一节中介绍的 generate() 函数进行配置。

运算符

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() 首先计算速度,然后限制速度值在 sminsmax 之间,最后计算新粒子位置。

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

操作符及其参数在工具箱中注册。粒子初始值在 [-100, 100] 范围内(pminpmax),并且速度在整个进化过程中限制在 [-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)

算法

一旦操作符在工具箱中注册,我们就可以通过首先创建一个新种群,然后应用原始的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基本示例可以在这里找到:examples/pso/basic

这是一段展示算法运行过程的视频,使用 matplotlib 绘制。红点代表迄今为止找到的最佳解决方案。

参考文献

[Poli2007]

Ricardo Poli, James Kennedy 和 Tim Blackwell, “粒子群优化概述”。群体智能。2007; 1: 33–57

[Knoek2003]

Arthur J. Knoek van Soest 和 L. J. R. Richard Casius, “并行遗传算法在解决复杂优化问题中的优势”。《生物力学工程杂志》。2003年; 125: 141–146