最优化算法之粒子群算法PSO

基本概念

粒子群算法(PSO)属于群智能算法的一种，1995 年由Eberhart 博士和kennedy 博士提出，源于对鸟群捕食的行为研究。该算法是通过模拟鸟群捕食行为设计的。假设区域里就只有一块食物（即通常优化问题中所讲的最优解），鸟群的任务是找到这个食物源。鸟群在整个搜寻的过程中，通过相互传递各自的信息，让其他的鸟知道自己的位置，通过这样的协作，来判断自己找到的是不是最优解，同时也将最优解的信息传递给整个鸟群，最终，整个鸟群都能聚集在食物源周围，即我们所说的找到了最优解，即问题收敛。
粒子群算法具有收敛速度快、参数少、算法简单易实现的优点（对高维度优化问题，比遗传算法更快收敛于最优解），但是也会存在陷入局部最优解的问题。

基本原理

位置和速度

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。鸟被抽象为没有质量和体积的微粒(点)，并延伸到N维空间，粒子i在N维空间的位置表示为矢量Xi＝(x1，x2，…，xN)，飞行速度表示为矢量Vi＝(v1，v2，…，vN)。每个粒子都有一个由目标函数决定的适应值(fitness value)，并且知道自己到目前为止发现的最好位置(pbest)和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置(gbest)(gbest是pbest中的最好值)，这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。

速度和位置的更新

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

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

公式(2)和 公式(3)被视为标准PSO算法。

算法流程

学习因子c1、c2分析

公式(2)和(3)中pbest和gbest分别表示微粒群的局部和全局最优位置。
当C1＝0时，则粒子没有了认知能力，变为只有社会的模型(social-only)：

称为全局PSO算法。粒子有扩展搜索空间的能力，具有较快的收敛速度，但由于缺少局部搜索，对于复杂问题
比标准PSO 更易陷入局部最优。
当C2＝0时，则粒子之间没有社会信息，模型变为只有认知(cognition-only)模型：

称为局部PSO算法。由于个体之间没有信息的交流，整个群体相当于多个粒子进行盲目的随机搜索，收敛速度慢，因而得到最优解的可能性小。

示例演示

用粒子群算法求解求解函数y=-x*(x-1) 在0,2上最大值。

#include 
#include 
#include 

using namespace std;

class AlgorithmPSO 
{
public:
	AlgorithmPSO()
	{

	}
	~AlgorithmPSO()
	{
		delete[] m_y;
		delete[] m_x;
		delete[] m_v;
		delete[] m_pbest;
	}

    //适应度计算函数，每个粒子都有它的适应度
     void fitnessFunction()
	 {
        for(int i=0;i<m_n;i++)
		{
            m_y[i] = -1.0f * m_x[i] * (m_x[i]-2.0);
        }
     }
    
	 //初始化
	 void init()
	 { 
        m_x = new double[m_n];
        m_v = new double[m_n];
        m_y = new double[m_n];
        m_pbest = new double[m_n];
        /***
         * 本来是应该随机产生的，为了方便演示，我这里手动随机落两个点，分别落在最大值两边
         */
        m_x[0] = 0.0;
        m_x[1] = 2.0;
        m_v[0] = 0.01;
        m_v[1] = 0.02;
        fitnessFunction();
        //初始化当前个体最优位置，并找到群体最优位置
        for(int i=0;i<m_n;i++)
		{
            m_pbest[i] = m_y[i];
            if(m_y[i]>gbest) gbest=m_y[i];
        }
       std::cout << ("Algorithm Starting , gbest:") << gbest << std::endl;
    }

    double getMAX(double a,double b)
	{
        return a>b?a:b;
    }

    //粒子群算法
    void PSO(int max)
	{
        for(int i=0;i<max;i++)
		{
            double w=0.4;
            for(int j=0;j<m_n;j++)
			{
                //更新位置和速度，下面就是我们之前重点讲解的两条公式。
				default_random_engine e;
				e.seed(time(0));
				std::normal_distribution<double> Normal_d(0, 1);
				double rd = Normal_d(e);

                m_v[j] = w * m_v[j]+c1 * rd * (m_pbest[j]-m_x[j]) + c2 * rd * (gbest-m_x[j]);
                if(m_v[j]>vmax) m_v[j] = vmax;//控制速度不超过最大值
                m_x[j] += m_v[j];
                
                //越界判断，范围限定在[0, 2]
                if(m_x[j]>2) m_x[j] = 2;
                if(m_x[j]<0) m_x[j] = 0;
                
            }
            fitnessFunction();
            //更新个体极值和群体极值
            for(int j=0;j<m_n;j++)
			{
                m_pbest[j]  =getMAX(m_y[j],m_pbest[j]);
                if(m_pbest[j]>gbest) gbest=m_pbest[j];
                std::cout << "Particle n" << j << ": x = " << m_x[j] <<"  " <<"v = " <<m_v[j] << std::endl;
            }
           std::cout  << (i+1)  << "iter , gbest = "  <<  gbest << std::endl;
        }
        
    }
private:
	int m_n = 2; //粒子个数，这里为了方便演示，我们只取两个，观察其运动方向
	double* m_y;
	double* m_x;
	double* m_v;
	double c1 = 2;
	double c2 = 2;
	double* m_pbest;
	double gbest;
	double vmax = 0.1; //速度最大值
};

int main()
{
	AlgorithmPSO ts;
	ts.init();
	ts.PSO(20);//为了方便演示，我们暂时迭代20次。
	return EXIT_SUCCESS;
}
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