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

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年提出];
  除了上述几种常见的群体智能算法以外,还有一些并不是广泛应用的群体智能算法,比如萤火虫算法、布谷鸟算法、蝙蝠算法以及磷虾群算法等等。

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

举个通俗的例子:

一群鸟在寻找食物,在这个区域内有一个食物,所有的鸟都不知道食物在哪里,但食物发出了香味,鸟通过香味的浓烈能判断出自己的当前位置距离食物有多远,同时鸟群之间是可以交流并告知离食物最近的鸟的位置的。试想一下这个时候会发生什么?

鸟甲:哈哈哈,原来我离食物最近!

鸟乙、丙、丁…:我得赶紧往鸟甲那里过去看看!

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

鸟某某:我刚才的位置好像离食物很近了,我得往那里再靠近点!

通过这样一个过程,不断的变化,最终鸟群会向食物方向聚集,达到目标。在这个变化的过程中,影响鸟的运动状态变化的有两个因素:

  • 离食物最近的鸟的位置
  • 自己之前达到过的离食物最近的位置

而鸟的每次的位置变化,除了考虑以上的两个因素还有一个因素就是惯性。而对位置的变化,通过变化量(速度)来表示。

通过以上的例子,我们可以把鸟抽象为没有质量和体积的微粒(点),并延伸到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)中,前期变化较慢,取值较大,维持了算法的全局搜索能力;后期变化较快,极大的提高了算法的局部寻优能力,通过这种权重的变化函数,可以取得比固定不变的那种情况要好的效果。

posted @ 2017-08-22 22:58  寒潭  阅读(10086)  评论(2编辑  收藏  举报