智能优化算法——粒子群算法原理与仿真程序

xiaozongpeng

写在前面：
粒子群算法很古老了，资料也很多，大部分寻优算法思想类似，本文稍作整理，并在文末分别附上C++、Python、MATLAB程序，如果本文对你有所帮助，请点赞支持一下！
这是遗传算法的原理与程序：
一、粒子群算法的概念

粒子群优化算法(PSO：Particle swarm optimization) 是一种进化计算技术（evolutionary computation）。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解。
PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。
二、粒子群算法分析

1、基本思想

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。
每个粒子在搜索空间中单独的搜寻最优解，记为个体当前极值，并将个体极值与整个粒子群里的其他粒子共享，整个粒子群的当前全局最优解为最优的那个个体极值。
粒子群中的所有粒子根据个体自身当前极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。下面的动图很形象地展示了PSO算法的过程：



PSO寻优算法过程示意

2、更新规则

PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。

在公式(1)、(2)中，i=1,2，。..，N，N是此粒子群中粒子总数。
：个体最优粒子
：全局最优粒子
  ：是粒子的速度
：介于（0,1）之间的随机数
：粒子的当前位置
和 ：学习因子，通常
的最大值为(大于0），如果大于  ，则 。
公式（1）、（2）为PSO的标准形式。

公式(1)的第一部分称为【记忆项】，表示上次速度大小和方向的影响；公式(1)的第二部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；公式(1)的第三部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。以上面两个公式为基础，形成了PSO的标准形式。

式中 叫做惯性因子，非负值。其值较大，全局寻优能力强，局部寻优能力弱；其值较小，全局寻优能力弱，局部寻优能力强。
动态  能获得比固定值更好的寻优结果。动态可在PSO搜索过程中线性变化，也可以根据PSO性能的某个测度函数动态改变。

目前采用较多的是线性递减权值（Linearly Decreasing Weight,LDW）策略。即：

：最大迭代次数
：初始惯性权值
：迭代至最大进化代数时的惯性权值
典型权值：

的引入，使PSO算法性能有了很大的提高，针对不同的搜索问题，可以调整全局和局部搜索能力，也使PSO算法能广泛地应用于很多实际问题。
公式(2)和 公式(3)被视为标准PSO算法。
3、PSO算法的流程和伪代码

粒子群算法是一种群体智能算法，通过跟随当前搜索到的最优值来寻找全局最优。实现容易、精度高、收敛快，在解决实际问题中具有很大的优越性。主要步骤可描述如下：

1、初始化粒子群位置和速度。
2、计算每个粒子的适应度，确定全局最优粒子gbest和个体最优粒子pbest。
3、判断算法收敛准则是否满足，若满足，则输出搜索结果；否则执行以下步骤。
4、根据gbest和pbest更新速度分量。
5、根据更新的速度分量更新粒子位置。
6、更新全局最优粒子和个体最优粒子。
7、返回步骤③。

4、PSO算法举例





三、粒子群算法(PSO)仿真程序

1、C++代码

#include <iostream>
#include <vector>
#include <cmath>
#include <map>
#include <algorithm>
#include <random>
#include <ctime>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

const int dim = 1;//维数
const int p_num = 10;//粒子数量
const int iters = 100;//迭代次数
const int inf = 999999;//极大值
const double pi = 3.1415;
//定义粒子的位置和速度的范围
const double v_max = 4;
const double v_min = -2;
const double pos_max = 2;
const double pos_min = -1;
//定义位置向量和速度向量
vector<double> pos;
vector<double> spd;
//定义粒子的历史最优位置和全局最优位置
vector<double> p_best;
double g_best;
//使用eigen库定义函数值矩阵和位置矩阵
Matrix<double, iters, p_num> f_test;
Matrix<double, iters, p_num> pos_mat;

//定义适应度函数
double fun_test(double x)
{
    double res = x * x + 1;
    return res;
}

//初始化粒子群的位置和速度
void init()
{
    //矩阵中所有元素初始化为极大值
    f_test.fill(inf);
    pos_mat.fill(inf);
    //生成范围随机数
    static std::mt19937 rng;
    static std::uniform_real_distribution<double> distribution1(-1, 2);
    static std::uniform_real_distribution<double> distribution2(-2, 4);
    for (int i = 0; i < p_num; ++i)
    {
        pos.push_back(distribution1(rng));
        spd.push_back(distribution2(rng));
    }
    vector<double> vec;
    for (int i = 0; i < p_num; ++i)
    {
        auto temp = fun_test(pos);//计算函数值
        //初始化函数值矩阵和位置矩阵
        f_test(0, i) = temp;
        pos_mat(0, i) = pos;
        p_best.push_back(pos);//初始化粒子的历史最优位置
    }
    std::ptrdiff_t minRow, minCol;
    f_test.row(0).minCoeff(&minRow, &minCol);//返回函数值矩阵第一行中极小值对应的位置
    g_best = pos_mat(minRow, minCol);//初始化全局最优位置
}

void PSO()
{
    static std::mt19937 rng;
    static std::uniform_real_distribution<double> distribution(0, 1);
    for (int step = 1; step < iters; ++step)
    {
        for (int i = 0; i < p_num; ++i)
        {
            //更新速度向量和位置向量
            spd = 0.5 * spd + 2 * distribution(rng) * (p_best - pos) +
                2 * distribution(rng) * (g_best - pos);
            pos = pos + spd;
            //如果越界则取边界值
            if (spd < -2 || spd > 4)
                spd = 4;
            if (pos < -1 || pos > 2)
                pos = -1;
            //更新位置矩阵
            pos_mat(step, i) = pos;
        }
        //更新函数值矩阵
        for (int i = 0; i < p_num; ++i)
        {
            auto temp = fun_test(pos);
            f_test(step, i) = temp;
        }
        for (int i = 0; i < p_num; ++i)
        {
            MatrixXd temp_test;
            temp_test = f_test.col(i);//取函数值矩阵的每一列
            std::ptrdiff_t minRow, minCol;
            temp_test.minCoeff(&minRow, &minCol);//获取每一列的极小值对应的位置
            p_best = pos_mat(minRow, i);//获取每一列的极小值，即每个粒子的历史最优位置
        }
        g_best = *min_element(p_best.begin(), p_best.end());//获取全局最优位置
    }
    cout << fun_test(g_best);
}

int main()
{
    init();
    PSO();
    system("pause");
    return 0;
}
2、Python代码

# -*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import math

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

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

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

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

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

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

def func(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(sizepop):
    pop = np.zeros((sizepop,2))
    v = np.zeros((sizepop,2))
    fitness = np.zeros(sizepop)

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

def getinitbest(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  

w = getweight()
lr = getlearningrate()
maxgen = getmaxgen()
sizepop = getsizepop()
rangepop = getrangepop()
rangespeed = getrangespeed()

pop,v,fitness = initpopvfit(sizepop)
gbestpop,gbestfitness,pbestpop,pbestfitness = getinitbest(fitness,pop)

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

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

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

        for j in range(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 = gbestfitness


plt.plot(result)
plt.show()3、MATLAB代码

(为缩减篇幅，这里给出主函数的m文件，如需全部matlab代码请评论区留言！)
clc
clear

%设置问题相关变量
popsize=2000;%个体数量
degree=100;%循环次数

[popmax,vmax]=Init(popsize);%初始化求最大值群体的位置和速度
[popmin,vmin]=Init(popsize);%初始化求最小值群体的位置和速度

%计算初始群体的适应度,群体最优个体
[value_max,gbestvaluemax,gbestmax]=calfitvaluemax(popmax);
[value_min,gbestvaluemin,gbestmin]=calfitvaluemin(popmin);

%将个体历史最优先设为初始化值
pbestmax=popmax;
pbestmin=popmin;

%创建图
set(gcf,'doublebuffer','on')
set(gcf,'Name','粒子群算法PSO演示')
axis([1 6 2 9])
grid on

for i=1:degree
    %产生新的位置
    [newpopmax,newvmax]=updatepop(popmax,vmax,pbestmax,gbestmax);
    [newpopmin,newvmin]=updatepop(popmin,vmin,pbestmin,gbestmin);
   
    %更新速度值
    vmax=newvmax;
    vmin=newvmin;
   
    %更新位置值
    popmax=newpopmax;
    popmin=newpopmin;
   
    %计算新的适应度值，该次循环中的群体最优
    [newvalue_max,newgbestvaluemax,newgbestmax]=calfitvaluemax(newpopmax);
    [newvalue_min,newgbestvaluemin,newgbestmin]=calfitvaluemin(newpopmin);
   
    %更新个体历史最优
    for j=1:popsize
        if newvalue_max(j)>value_max(j)
            pbestmax(:,j)=newpopmax(:,j);
        end
        if newvalue_min(j)<value_min(j)
            pbestmin(:,j)=newpopmin(:,j);
        end      
    end
   
    %更新群体最优
    if newgbestvaluemax>gbestvaluemax
        gbestvaluemax=newgbestvaluemax;
        gbestmax=newgbestmax;
    end
   
    if newgbestvaluemin<gbestvaluemin
        gbestvaluemin=newgbestvaluemin;
        gbestmin=newgbestmin;
    end
   
    %绘制动态运动散点图
    plot(newpopmax(1,:),newpopmax(2,:),'r*',newpopmin(1,:),newpopmin(2,:),'go');
    drawnow
   
end
四、参考资料

1.最优化算法之粒子群算法（PSO）_青萍之末的博客-CSDN博客_粒子群算法
2.EddyGao - Overview

如果本文对你有所帮助，请点赞支持一下！

Baste

想要matlab代码