粒子群优化算法

资料来源

粒子群优化算法_百度百科 (baidu.com)

简介

定义

粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法、或微粒群优化算法。是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。它可以被纳入多主体优化系统(Multiagent Optimization System, MAOS)。

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

模拟捕食

PSO模拟鸟群的捕食行为。一群鸟在随机搜索食物，在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻离食物最近的鸟的周围区域。

启示

PSO从这种模型中得到启示并用于解决优化问题。PSO中，每个优化问题的解都是搜索空间中的一只鸟。我们称之为“粒子”。所有的粒子都有一个由被优化的函数决定的适应值(fitnessvalue)，每个粒子还有一个速度决定他们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。

PSO初始化

PSO初始化为一群随机粒子(随机解)，然后通过迭代找到最优解，在每一次迭代中，粒子通过跟踪两个“极值”来更新自己。第一个就是粒子本身所找到的最优解，这个解叫做个体极值pBest，另一个极值是整个种群找到的最优解，这个极值是全局极值gBest。另外也可以不用整个种群而只是用其中一部分最优粒子的邻居，那么在所有邻居中的极值就是局部极值。

基本思想

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

算法介绍

在找到这两个最优值时,粒子根据如下的公式来更新自己的速度和新的位置

img

(a)

img

(b)

v[] 是粒子的速度，present[] 是当前粒子的位置。pbest[] 和 gbest[] 如前定义。rand() 是介于(0,1)之间的随机数。c1，c2是学习因子。通常c1=c2=2。

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

另一种解释：

这里写图片描述

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

这里写图片描述

img

伪代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
For each particle

____Initialize particle

END

Do

____For each particle

________Calculate fitness value

________If the fitness value is better than the best fitness value (pBest) in history

____________set current value as the new pBest

____End

____Choose the particle with the best fitness value of all the particles as the gBest

____For each particle

________Calculate particle velocity according equation (a)

________Update particle position according equation (b)

____End
While maximum iterations or minimum error criteria is not attained
在每一维粒子的速度都会被限制在一个最大速度Vmax，如果某一维更新后的速度超过用户设定的Vmax，那么这一维的速度就被限定为Vmax。

PSO算法案例

这里写图片描述

这里写图片描述

代码demo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#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[i]);//计算函数值
        //初始化函数值矩阵和位置矩阵
        f_test(0, i) = temp;
        pos_mat(0, i) = pos[i];
        p_best.push_back(pos[i]);//初始化粒子的历史最优位置
    }
    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[i] = 0.5 * spd[i] + 2 * distribution(rng) * (p_best[i] - pos[i]) +
                2 * distribution(rng) * (g_best - pos[i]);
            pos[i] = pos[i] + spd[i];
            //如果越界则取边界值
            if (spd[i] < -2 || spd[i] > 4)
                spd[i] = 4;
            if (pos[i] < -1 || pos[i] > 2)
                pos[i] = -1;
            //更新位置矩阵
            pos_mat(step, i) = pos[i];
        }
        //更新函数值矩阵
        for (int i = 0; i < p_num; ++i)
        {
            auto temp = fun_test(pos[i]);
            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[i] = 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;
}