首页 > 模式识别 > 粒子群优化(Particle Swarm Optimization, PSO)

粒子群优化(Particle Swarm Optimization, PSO)

2013年7月16日 发表评论 阅读评论

今天突然开始发闲了，因为老板好像出差了，然后今天下午原计划是去另一个校区搞一些比较重要的事情的，但是我问那边老师下午在么？他在1点50分才回我短信说。。。2点之前他会在的。。。然后下午突然就觉得好空虚，没事干了。。

好吧，所以，今天继续上来补一篇比较简单的算法文吧，应该没什么人会看不懂的吧。。。估计。。。

RT所示，叫做粒子群优化,因为我一直管叫PSO，中文名在我脑海中的印象一直是“粒子XX算法”，我今天才突然意识到，原来有优化这两个字的，然后我就在想，这个叫做“优化”的算法都那么简单了，那么那个叫做“粒子算法”的东西是有多简单啊。。。然后查了一下，好吧，没有“粒子算法”这么一说，或者说，“粒子算法”就是“粒子群优化算法”本身。。。

可能有人听过蚁群算法，觉得粒子群算法会不会跟蚁群算法差不多，其实。。不是差得多不多的问题，完全就是两个东西，生态一点的说法，蚁群算法是模拟蚂蚁找回家的路这种最优路径问题，但是粒子群模拟的是鸟类找最适合栖息地的过程，虽说都是找最优解的问题，但是你敢说鸟类的做法可以放到蚂蚁身上？算了，过阵子我有空就再搞一篇蚁群算法的吧。。。这里姑且先开这个坑。。。好吧，虽说是模仿鸟类，但是为了显得“专业”，我们还是管鸟类叫做“粒子”好了，鸟群叫做“粒子群”！

既然用了粒子来表示鸟类，那么粒子就要有一个基本特性！要会飞！！怎么描述会飞这件事情，神说，要有个速度，所以在粒子群算法中，我们赋予每个粒子一个速度(高中物理告诉我们，速度是矢量！！)。

既然有了速度，就要有位置这个概念，速度才得以成立，所以，每个粒子每一时刻都有一个位置，每一时刻的位置取决于前一时刻的位置和速度。

再然后，既然是最优化问题，那么我们应该赋予每一个位置一个权值，但是为了显得专业(因为很重要，所以要说两遍！！)，那个叫适应度(这个词在大部分寻找最优值的算法里面都会被用上，而估计起源是遗传算法。。没考证过。。)。

然后既然说了这个算法是粒子群协同合作找最优值(也就是找适应度最大的点)，那么就要有信息的共享，为了简单说明，我们可以像魔兽争霸或者DOTA(重要声明：我是不打DOTA的！！争霸一生推！！)一般，虽然是友军的视野，我们也可以看到，也就是说，每个粒子的信息(比如它的位置，速度，最重要的是当前的适应度)都是被所有粒子所共享的，啊，这么说来，争霸/DOTA这个例子还举得不好，应该举御坂网络才对的。。。

虽说是信息共享，但是对于这个算法中有用的只有一个，所有粒子都知道目前所有粒子所找到的最优的那个点在哪里，你们可以想象成每当某个粒子找到某个位置比目前找到的最好的还要好的时候，它就会大吼一句“德玛西亚！！”，然后把坐标的位置告诉全部的小伙伴。（之前自以为自己找到的是最优点的粒子和它的小伙伴都惊呆了！！）

然后最后一个，就是每个粒子都要刷一个人生成就，就是自己找到的目前的最优位置，这个它会时时刻刻死死地记住的。

总之，要说明的就这些，然后描述一下算法：

描述之前，老师告诉我们，先要定义变量，首先是目前为止找到的全局最优点的位置Pg和最优点的适应度Fg，然后是每个粒子自己找到的最优点的位置Pi和最优点的适应度Fi，然后每个粒子的当前速度Vi,每个粒子当前的位置Si，额，好像就没了。。【所以说这个算法很简单啊。。】

  1. 随机初始化所有粒子的初始位置，速度，计算出当前位置的适应度，作为每个粒子的Pi和Fi,然后找出最大的Fi作为Fg，并且这个最大的Fi的位置作为Pg。
  2. 更新每个粒子的速度和位置。(更新方法，显然的。。在后文。。)
  3. 计算新的的Pi和Fi
  4. 根据步骤1的方法，更新Fg和Pg
  5. 判断是否达到迭代停止条件(比如达到迭代次数上限)，如果没达到，回到2，不然，就到6
  6. 结束，该干嘛干嘛去。。
然后就是给出第二步更新算法，公式如下：
 Vi_new=ω Vi_old+C1×rand()×(Pi-Si)+C2×rand()×(Pg-Si)
Si_new= Si_old+Vi_new
其中要说明的三个参数，ω，C1和C2，都是常数，据论文和维基娘所述，ω叫惯性权重(inertia weight)，而C1和C2，叫加速常数。然后rand()生成的随机数，一般是[0,1]这个范围

公式很容易理解，就是每个粒子即会一定程度上因为惯性保持现有的速度，觉得自己要保持自我，为群体找到更好的栖息地，又一定程度上会往自己找到的最优点那里飞，想着我之前那个最好的栖息地附近会不会有更好的点呢？又会一定程度上往群体找到的最好的点那里偏移，因为，你懂得嘛，随大流~最后因为这种极度矛盾的心理挣扎之后，出现了公式所述的速度更新公式。。。【好吧，一般装B的论文里面都是管第二项叫自我认知，第三项叫社会认知的。。。但是我觉得。。。这种说法好中二的感觉。。】

第二个公式。。。应该不用解释吧。。。看不懂就在Vi_new后面乘个Δt，但是因为为了离散化处理，所以Δt都是1而已。。。

以上，就是正篇的全文了。。。

然后这里补充点东西，但是不做过多说明，属于。。。超纲的。。

一个就是，我刚刚就很想吐槽DOTA里面共享友军视野这种事情，凭什么你队友看到的东西你要和自己亲眼看到一样。。虽然游戏里面不这么做，必然是不行的，但是最为仿生的算法，就算你的一句“德玛西亚”在大声，也不可能传遍整个大江南北的，你以为你在中国好声音现场啊，就算你在中国好声音现场，你以为每个粒子都有一段悲惨的过去啊，中国好故事不是什么人都有的。。所以就出现了一种局部共享信息的算法，就是说，对于每个粒子的Pg都是不一样的，只是自己周围能“听到”的朋友告诉的最优的位置。这种做法的好处是什么呢？其实很容易想明白，就是全局信息共享的话，很容易陷入局部最优点，也就是全部粒子都一窝蜂的往全局最优点那里飞，局部共享的话可以很好的克服这个问题。但是，相对熟悉各种人工智能算法分析的童鞋都知道，只要一个算法通过改进，使得可以可以很好的避开局部最优点，那么付出的代价里面必然是。。。收敛速度。。

另一个想说的就是，以前在搞数模的时候看到的，其实粒子群可以用来解决01组合这种问题，就是把适应度和位置全部量化成0和1，只有速度还是连续变化的，具体的请去找相关专业文献看看，因为我也没写过这个算法。。

然后下面是我以前写的一个模拟二维粒子群找山峰的最高值的演示程序(代码在文中最后)，背景里面越白的地方代表那个值越大，然后看看找的过程(注意，其实视屏是可以调超清的。。优酷有没有广告我就不知道了)。。。

为什么突然想放这个上来呢？因为。。我觉得视频最后那一段很壮观的样子诶。。。

下面是代码，很不好意思的一点是。。。以前写的，显示时用OPENCV做的。。。所以没装的电脑都跑不了。。。所以看看视频就好了~

Code:

#ifdef _DEBUG
#pragma comment ( lib, "cxcore200d.lib" )
#pragma comment ( lib, "cv200d.lib" )
#pragma comment ( lib, "highgui200d.lib" )
#else
#pragma comment ( lib, "cxcore200.lib" )
#pragma comment ( lib, "cv200.lib" )
#pragma comment ( lib, "highgui200.lib" )
#endif
#include "cv.h"
#include "highgui.h"

#include <iostream>
#include <stdlib .h>
#include </stdlib><stdlib .h>
#include <time .h>
using namespace std;

#define WINDOW_NAME "PSO"
#define MAX_PARTICAL_NUM 100
#define MAX_ITERATION 10000
#define FIELD_HEIGTH 500
#define FIELD_WIDTH 600
#define DT 1
#define CENTER_NUM 5

class str_vec
{
public:
   str_vec(){};
   str_vec(double xx,double yy):x(xx),y(yy){};
   double x;
   double y;
};

struct str_partical
{
   str_vec pos;
   str_vec velocity;
   double fitness;

   str_vec best_pos;
   double best_fitness;
};

str_partical particle[MAX_PARTICAL_NUM];
str_partical best_particle;

IplImage* pic_BackGround = NULL;
IplImage* frame = NULL;

int** Fitness_board;
const int CENTER_X[CENTER_NUM]={
                               FIELD_WIDTH/4,FIELD_WIDTH/4,FIELD_WIDTH*3/4,FIELD_WIDTH*3/4,FIELD_WIDTH/2};
const int CENTER_Y[CENTER_NUM]={
                               FIELD_HEIGTH/4,FIELD_HEIGTH*3/4,FIELD_HEIGTH/4,FIELD_HEIGTH*3/4,FIELD_HEIGTH/2};
double weigth[CENTER_NUM] = {2.1,2.4,2.3,2,2.9};
static void InitFitness()
{
   Fitness_board = new int*[FIELD_HEIGTH];
   for(int i = 0;i < FIELD_HEIGTH;i++)
   {
       Fitness_board[i] = new int[FIELD_WIDTH];
       for(int j = 0;j < FIELD_WIDTH;j++)
       {
           Fitness_board[i][j] = 0;
       }
   }

   double maxx = 0;
   for(int i = 0;i < FIELD_HEIGTH;i++)
   {
       for(int j = 0;j < FIELD_WIDTH;j++)
       {
           for(int k = 0;k < CENTER_NUM;k++)
           {
               double temp1 = double((CENTER_X[k]-j)*(CENTER_X[k]-j)+(CENTER_Y[k]-i)*(CENTER_Y[k]-i))/3600;
               Fitness_board[i][j] += (weigth[k]*1024/(1+temp1));
           }
           if(Fitness_board[i][j] > maxx)
           {
               maxx = Fitness_board[i][j];
           }
       }
   }
   for(int i = 0;i < FIELD_HEIGTH;i++)
   {
       for(int j = 0;j < FIELD_WIDTH;j++)
       {
           Fitness_board[i][j] *= (255/maxx);
       }
   }
}

static double GetFitness(str_vec pos)
{
   if(pos.x <= 0 || pos.x >= FIELD_WIDTH || pos.y < = 0 || pos.y >= FIELD_HEIGTH)
       return 0;
   return Fitness_board[int(pos.y)][int(pos.x)];
}

static void InitParticle(int i)
{
   particle[i].pos.x = rand()%FIELD_WIDTH;
   particle[i].pos.y = rand()%FIELD_HEIGTH;
   particle[i].velocity.x = (rand()%500-500)/100;
   particle[i].velocity.y = (rand()%500-500)/100;
   particle[i].fitness = GetFitness(particle[i].pos);

   particle[i].best_pos = particle[i].pos;
   particle[i].best_fitness = particle[i].fitness;
}

static void InitAllParticle()
{
   InitParticle(0);
   best_particle = particle[0];
   for(int i = 1;i < MAX_PARTICAL_NUM;i++)
   {
       InitParticle(i);
       if(particle[i].fitness > best_particle.fitness)
       {
           best_particle = particle[i];
       }
   }
}

static void InitBackGround()
{
   pic_BackGround = cvCreateImage(cvSize(FIELD_WIDTH,FIELD_HEIGTH),IPL_DEPTH_8U,3);
   frame = cvCreateImage(cvSize(FIELD_WIDTH,FIELD_HEIGTH),IPL_DEPTH_8U,3);
   for(int i = 0;i < FIELD_HEIGTH;i++)
   {
       for(int j = 0;j < FIELD_WIDTH;j++)
       {
           CV_IMAGE_ELEM(pic_BackGround,uchar,i,3*j) =
               CV_IMAGE_ELEM(pic_BackGround,uchar,i,3*j+1) =
                   CV_IMAGE_ELEM(pic_BackGround,uchar,i,3*j+2) = uchar(GetFitness(str_vec(j,i)));
       }
   }
}

static void UpdateParticle(int i)
{
   particle[i].pos.x += particle[i].velocity.x * DT;
   particle[i].pos.y += particle[i].velocity.y * DT;

   particle[i].fitness = GetFitness(particle[i].pos);

//     if(abs(best_particle.best_pos.y - particle[i].pos.y) + abs(best_particle.best_pos.x - particle[i].pos.x) < 20)
//     {
//         particle[i].velocity.x /= 1.01;
//         particle[i].velocity.y /= 1.01;
//     }
//     else
   {
       particle[i].velocity.x += (
           double(rand()%100)/99999*(particle[i].best_pos.x - particle[i].pos.x) +
           double(rand()%100)/99999*(best_particle.best_pos.x - particle[i].pos.x)
           );
       particle[i].velocity.y += (
           double(rand()%100)/99999*(particle[i].best_pos.y - particle[i].pos.y) +
           double(rand()%100)/99999*(best_particle.best_pos.y - particle[i].pos.y)
           );
   }
   if (particle[i].fitness >= particle[i].best_fitness)
   {
       particle[i].best_fitness = particle[i].fitness;
       particle[i].best_pos = particle[i].pos;
   }

   if (particle[i].best_fitness >= best_particle.best_fitness)
   {
       cvLine(pic_BackGround,cvPoint(best_particle.best_pos.x,best_particle.best_pos.y),
           cvPoint(particle[i].best_pos.x,particle[i].best_pos.y),cvScalar(255,0,255));
       best_particle.best_pos = particle[i].best_pos;
       best_particle.best_fitness = particle[i].best_fitness;
   }
}

static void UpdatePSO()
{
   for(int i = 0;i < MAX_PARTICAL_NUM;i++)
   {
       UpdateParticle(i);
   }
}

static void DrawParticle()
{
   cvCopyImage(pic_BackGround,frame);
   for(int i = 0;i < MAX_PARTICAL_NUM;i++)
   {
       if( particle[i].pos.x <= 0 ||
           particle[i].pos.x >= FIELD_WIDTH ||
           particle[i].pos.y < = 0 ||
           particle[i].pos.y >= FIELD_HEIGTH)
           continue;

       cvCircle(frame,cvPoint(particle[i].pos.x,particle[i].pos.y),1,cvScalar(0,0,255));
       cvLine(frame,cvPoint(particle[i].pos.x,particle[i].pos.y),
               cvPoint(particle[i].pos.x+particle[i].velocity.x,particle[i].pos.y+particle[i].velocity.y),
               cvScalar(255,0,0));
   }
   cvCircle(frame,cvPoint(best_particle.best_pos.x,best_particle.best_pos.y),5,cvScalar(0,255,0));
}

int main(int argc,char* argv[])
{
   srand((unsigned int)time(0));
   InitFitness();
   InitAllParticle();
   InitBackGround();
   DrawParticle();
   cvNamedWindow(WINDOW_NAME);
   cvShowImage(WINDOW_NAME,frame);
   cvWaitKey(0);

   for(int i = 0;i < MAX_ITERATION;i++)
   {
       UpdatePSO();
       DrawParticle();
       cvShowImage(WINDOW_NAME,frame);
       cvWaitKey(20);
   }
   cout<<"pos:"<<best_particle.best_pos.x<<" "<<best_particle.best_pos.y<<endl;
   cout<<"fitness:"<<best_particle.best_fitness<<endl;
   cvWaitKey(0);
   return 0;
}

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com

分类: 模式识别 标签: , , ,
  1. 静水流深
    2015年12月23日16:33 | #1
    回复 | 引用

    哈哈，我在知网还有一篇水文章是关于粒子群优化的

验证码：9 + 8 = ?

友情提示：留言可以使用大部分html标签和属性;

添加代码示例:[code lang="cpp"]your code...[/code]

添加公式请用Latex代码，前后分别添加两个$$