【智能算法】粒子群寻优算法

时间  2019-12-09
标签 智能算法 粒子 算法 繁體版
原文   原文链接

1.理论基础

粒子群算法（particle swarm optimization，PSO）是计算智能领域中的一种生物启发式方法，属于群体智能优化算法的一种，常见的群体智能优化算法主要有以下几类：
　　（1）蚁群算法（Ant Colony Optimization，简称ACO）[1992年提出]；
　　（2）粒子群优化算法（Particle Swarm Optimization，简称PSO）[1995年提出]（简单易于实现，也是目前应用最为普遍的群体智能优化算法）；
　　（3）菌群优化算法（Bacterial Foraging Optimization，简称BFO）[2002年提出]；
　　（4）蛙跳算法（Shuffled Frog Leading Algorithm，简称SFLA）[2003年提出]；
　　（5）人工蜂群算法（Artificial Bee Colony Algorithm，简称ABC）[2005年提出]；
　　除了上述几种常见的群体智能算法之外，还有一些并非普遍应用的群体智能算法，好比萤火虫算法、布谷鸟算法、蝙蝠算法以及磷虾群算法等等。python

而其中的粒子群优化算法（PSO）源于对鸟类捕食行为的研究，鸟类捕食时，找到食物最简单有限的策略就是搜寻当前距离食物最近的鸟的周围。算法

举个通俗的例子：dom

一群鸟在寻找食物，在这个区域内有一个食物，全部的鸟都不知道食物在哪里，但食物发出了香味，鸟经过香味的浓烈能判断出本身的当前位置距离食物有多远，同时鸟群之间是能够交流并告知离食物最近的鸟的位置的。试想一下这个时候会发生什么？ide

鸟甲：哈哈哈，原来我离食物最近！函数

鸟乙、丙、丁…:我得赶忙往鸟甲那里过去看看！学习

同时各自鸟在位置不停变化的时候，离食物的距离也不断的变化，因此必定有过离食物最近的位置，这是他们的一个参考。优化

鸟某某：我刚才的位置好像离食物很近了，我得往那里再靠近点！spa

经过这样一个过程，不断的变化，最终鸟群会向食物方向汇集，达到目标。在这个变化的过程当中，影响鸟的运动状态变化的有两个因素：3d

  • 离食物最近的鸟的位置
  • 本身以前达到过的离食物最近的位置

而鸟的每次的位置变化，除了考虑以上的两个因素还有一个因素就是惯性。而对位置的变化，经过变化量（速度）来表示。code

经过以上的例子，咱们能够把鸟抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子I 在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)．每一个粒子都有一个由目标函数决定的适应值(fitness value)，而且知道本身到目前为止发现的最好位置(pbest)和如今的位置Xi．这个能够看做是粒子本身的飞行经验．除此以外，每一个粒子还知道到目前为止整个群体中全部粒子发现的最好位置(gbest)(gbest是pbest中的最好值)．这个能够看做是粒子同伴的经验．粒子就是经过本身的经验和同伴中最好的经验来决定下一步的运动。归结为下面的两个公式：

                                                    (公式-1）

                                                                                              （公式-2）

其中，为惯性权重因子，（什么是权重因子，这个后面会讲到，此处大概知道有这么一个参数便可）；d=1,2，…，D；i=1,2，…，n；k为当前迭代次数；Vid为粒子的速度；c1和c2是非负的常数，称为加速因子；r1和r2是分布于[0,1]区间的随机数。为了防止粒子的满目搜索，通常建议限制位置和速度在必定的区间。

2.粒子群寻优算法的案例与实现

假设有以下非线性函数，须要在必定范围内找出最大值和最大值的位置，本文用Python实现。

                                                    （公式-3）

咱们先经过python把上面函数使用matplotlib把它画处理，让咱们有个直观的感觉。

# -*- coding: utf-8 -*-
from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(-2, 2, 0.05)
Y = np.arange(-2, 2, 0.05)
X, Y = np.meshgrid(X, Y)
Z = np.sin(np.sqrt(X**2+Y**2))/np.sqrt(X**2+Y**2)+np.exp((np.cos(2*np.pi*X)+np.cos(2*np.pi*Y))/2)-2.71289

ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='hot')
plt.show()

此处使用mpl_toolkits.mplot3d画3D图，在代码中，分别使用numpy中的arange建立等差长度为0.05的数据点向量，区间为[-2,2]之间，而后使用meshgrid将向量转换给对应的矩阵，这样的话，后面的函数式中就可使用np中的三角函数和数学函数了，若是不这样的话，还得一个一个点的根据公式去作循环，拼接Z的向量点。

运行出来的图以下：

figure_1

从函数图形中能够看出，这个函数有不少的局部极大值，而真正的极限位置就是在（0,0），极大值为1.005。下面就经过粒子群算法来找寻这个极大值和极大值的位置。根据粒子群算法函数求极值的算法流程以下图所示：

image

在pycham做为python的IDE中，创建一个名为PSO.py的文件，文件中编写一个名PSO的类。

2.1 参数的设定

设置惯性参数为常数1，也就是公式-1中的为常数1；设置加速度因子为非负的两个常数，公式-1的中的c1和c2；设置迭代次数为300；设置粒子群的规模为20；为了防止粒子的盲目搜索，将位置和速度限制在区间[-2,2]、[-0.5,0.5]；

具体代码以下：

def __init__(self):
        self.w = self.getweight()
        self.lr = self.getlearningrate()
        self.maxgen = self.getmaxgen()
        self.sizepop = self.getsizepop()
        self.rangepop = self.getrangepop()
        self.rangespeed = self.getrangespeed()

    def getweight(self):
        # 惯性权重
        weight = 1
        return weight

    def getlearningrate(self):
        # 分别是粒子的个体和社会的学习因子，也称为加速常数
        lr = (0.49445,1.49445)
        return lr

    def getmaxgen(self):
        # 最大迭代次数
        maxgen = 300
        return maxgen

    def getsizepop(self):
        # 种群规模
        sizepop = 20
        return sizepop

    def getrangepop(self):
        # 粒子的位置的范围限制,x、y方向的限制相同
        rangepop = (-2,2)
        return rangepop

    def getrangespeed(self):
        # 粒子的速度范围限制
        rangespeed = (-0.5,0.5)
        return rangespeed

2.2 种群初始化

随机初始化粒子位置和粒子的速度，并根据适应度函数计算粒子适应度值，咱们的适应度函数就是公式-3这个函数。

首先定义一个适应度函数

def func(self,x):
        # x输入粒子位置
        # y 粒子适应度值
        if (x[0]==0)&(x[1]==0):
            y = np.exp((np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1]))/2)-2.71289
        else:
            y = np.sin(np.sqrt(x[0]**2+x[1]**2))/np.sqrt(x[0]**2+x[1]**2)+np.exp((np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1]))/2)-2.71289
        return y

而后初始化粒子和计算适应度值

def initpopvfit(self,sizepop):
        pop = np.zeros((sizepop,2))
        v = np.zeros((sizepop,2))
        fitness = np.zeros(sizepop)

        for i in range(sizepop):
            pop[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]
            v[i] = [(np.random.rand()-0.5)*self.rangepop[0]*2,(np.random.rand()-0.5)*self.rangepop[1]*2]
            fitness[i] = self.func(pop[i])
        return pop,v,fitness

2.3 寻找初始化后的极值

def getinitbest(self,fitness,pop):
        # 群体最优的粒子位置及其适应度值
        gbestpop,gbestfitness = pop[fitness.argmax()].copy(),fitness.max()
        #个体最优的粒子位置及其适应度值,使用copy()使得对pop的改变不影响pbestpop，pbestfitness相似
        pbestpop,pbestfitness = pop.copy(),fitness.copy()

        return gbestpop,gbestfitness,pbestpop,pbestfitness

此处须要分别找到群体极值和个体粒子的极值。

2.4 迭代寻优

经过循环迭代，不断的更新粒子的位置和速度，根据新粒子的适应度值更新个体和群体的极值。这部分的代码以下：

def run(self):
        pop, v, fitness = self.initpopvfit(self.sizepop)
        gbestpop, gbestfitness, pbestpop, pbestfitness = self.getinitbest(fitness, pop)

        result = np.zeros(self.maxgen)
        for i in range(self.maxgen):
            # 速度更新
            for j in range(self.sizepop):
                v[j] += self.lr[0] * np.random.rand() * (pbestpop[j] - pop[j]) + self.lr[1] * np.random.rand() * (
                gbestpop - pop[j])
            v[v < self.rangespeed[0]] = self.rangespeed[0]
            v[v > self.rangespeed[1]] = self.rangespeed[1]

            # 粒子位置更新
            for j in range(self.sizepop):
                pop[j] += 0.5 * v[j]
            pop[pop < self.rangepop[0]] = self.rangepop[0]
            pop[pop > self.rangepop[1]] = self.rangepop[1]

            # 适应度更新
            for j in range(self.sizepop):
                fitness[j] = self.func(pop[j])

            for j in range(self.sizepop):
                if fitness[j] > pbestfitness[j]:
                    pbestfitness[j] = fitness[j]
                    pbestpop[j] = pop[j].copy()

            if pbestfitness.max() > gbestfitness:
                gbestfitness = pbestfitness.max()
                gbestpop = pop[pbestfitness.argmax()].copy()

            result[i] = gbestfitness
        return result

2.4 结果分析

经过在主函数文件中，实例化PSO类，调用run这个函数，运行上面的代码

from mpl_toolkits.mplot3d  import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from PSO import *

pso = PSO()
result = pso.run()
plt.figure()
plt.plot(result)
plt.show()

当迭代了300次以后，画出每代最优个体适应度值的曲线，结果以下：

figure_2

最终获得的最优个体适应度值为：1.00538866933，对应的粒子位置为：[  2.37513360e-05   3.41265823e-04]，经过比较，PSO算法能获得最优值，接近函数实际的最优值1.005。

经过以上，证实PSO算法具备较强的函数极值寻优能力。

3.延伸与算法改进

讲到这里，咱们再来回顾一下以前粒子更新速度（公式-1）中的这个惯性权重因子，所谓惯性权重因子，体现的是粒子继承上一次迭代速度的能力。较大的一个惯性权重因子有利于全局的搜索，而一个小的惯性权重因子则更加利于局部的搜索。为了更好的平衡算法的全局搜索与局部搜索的能力，常给惯性权重因子改成随迭代次数变化的一个函数，通常经常使用的惯性权重因子函数有如下几种：

          （公式-4）

4                                   （公式-5）

2                       （公式-6）

3                 （公式-7）

 

这几种的变化图以下：

image

从变化图中，能够看到（公式-5）中，前期变化较慢，取值较大，维持了算法的全局搜索能力；后期变化较快，极大的提升了算法的局部寻优能力，经过这种权重的变化函数，能够取得比固定不变的那种状况要好的效果。

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程