单目标 PSO（粒子群算法）

PSO也写了一个多月了吧，记得是今年3月底完成了

差不多是破釜沉舟式写出来的，当时写出来那个激动哇。。

因为就花一天，而且很简单，更主要的是，写过了JADE，这个写起来就思路清晰很多了。

我觉得这个问题是比较有意思的一个问题，它是模拟鸟群捕食行为的一种算法：

设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻目前离食物最近的鸟的周围区域。

 

但是这个问题又继而演化成了一个物理问题了（这里得打个问号了，是我自己这么认为的）

PSO最主要的是有一个速度公式和一个位移公式，如下：

v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[])

present[] = persent[] + v[]

这里面的3个参数w,c1,c2是非常之重要的，w是惯性权重，如果它太大，相当于运动的惯性系数太大，会突然间超过目标（停不下来。。。），如果太小，那么飞行到目标解又得花很长的时间。更奇葩的应该是c1,c2的选取吧，等下看我程序里给的值吧（其实这里，标准算法的值是w=1.0,c1=c2=2.0的）

 

先看下PSO的基本流程吧

再贴一下百度上找的伪代码：

clear all;

clc;

format long;

%------给定初始化条件----------------------------------------------

c1=1.4962;             %学习因子1

c2=1.4962;             %学习因子2

w=0.7298;              %惯性权重

MaxDT=1000;            %最大迭代次数

D=10;                  %搜索空间维数（未知数个数）

N=40;                  %初始化群体个体数目

eps=10^(-6);           %设置精度(在已知最小值时候用)

%------初始化种群的个体(可以在这里限定位置和速度的范围)------------

for i=1:N

    for j=1:D

        x(i,j)=randn;  %随机初始化位置

        v(i,j)=randn;  %随机初始化速度

    end

end

%------先计算各个粒子的适应度，并初始化Pi和Pg----------------------

for i=1:N

    p(i)=f(x(i,:),D);

    y(i,:)=x(i,:);

end

pg=x(1,:);             %Pg为全局最优

for i=2:N

    if f(x(i,:),D)<f(pg,D)

        pg=x(i,:);

    end

end

%------进入主要循环，按照公式依次迭代，直到满足精度要求------------

for t=1:MaxDT

    for i=1:N

        v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));

        x(i,:)=x(i,:)+v(i,:);

        if f(x(i,:),D)<p(i)

            p(i)=f(x(i,:),D);

            y(i,:)=x(i,:);

        end

        if p(i)<f(pg,D)

            pg=y(i,:);

        end

    end

    Pbest(t)=f(pg,D);

end

%------最后给出计算结果

disp('*************************************************************')

disp('函数的全局最优位置为：')

Solution=pg'

disp('最后得到的优化极值为：')

Result=f(pg,D)

disp('*************************************************************')

%------算法结束---DreamSun GL & HF-----------------------------------

最后再贴下我自己写的代码吧（测试函数依旧用的JADE里的第一个，精度也非常精确了~~）

代码风格依旧很挫。。。。

#include<stdio.h>
#include<time.h>
#include<stdlib.h>
#include<iostream>
//#include<random>
#include<cmath>
#include<malloc.h>
using namespace std;


#define URAND rand()/(RAND_MAX+1.0)


const int P_num=100;
const double PI=acos(-1.0);
const int G=1500;
const double w=0.7162; //用0.72//// 崔文明用的0.7162，因此我也改用了，就相当准确了，不过标准的是0.729
                       //可以考虑线性变化的w的取值...
                       //w=0.9-t/T*0.5,不过测试之后，效果更差了
const double c1=1.49445;
const double c2=1.49445;
const int D=30;
int times;
struct individual
{
    double p[D+10];
    double v[D+10];
    double fitness;
    double pbest[D+10];
    double val;
} gbest,P[P_num+10];

double rand_low_high(double low,double high)
{
    return (high-low)*URAND+low;
}
double calc_val(individual P)
{
    double val=0;
    for(int i=1;i<=D;i++)
    val+=P.p[i]*P.p[i];
    return val;
}
void calc_P_val()
{
    for(int i=1; i<=P_num; i++)
    {
        P[i].val=calc_val(P[i]);
    }
}
void init_P()//初始化
{

    srand((unsigned)time(NULL));
    for(int i=1; i<=P_num; i++)
    {
        for(int j=1; j<=D; j++)
        {
            P[i].p[j]=rand_low_high(-100,100);
            P[i].v[j]=rand_low_high(-100,100);
        }
    }

    calc_P_val();

    //找出第一代里的pbest和gbest
    for(int i=1;i<=P_num;i++)
    {
        for(int j=1;j<=D;j++)
        P[i].pbest[j]=P[i].p[j];
        P[i].fitness=P[i].val;
    }
    gbest=P[1];
    //for(int i=2;i<=P_num;i++)
    //if(P[i].fitness>gbest.fitness)
    //gbest=P[i];
}

int main()
{
    srand((unsigned)time(NULL));
    init_P();
 //   double w=0.7;
    for(times=1; times<=G; times++)//其实更新应该是可以从times=2开始的
    //for(times=2; times<=G; times++)
    {
        for(int j=1;j<=P_num;j++)
        {
            for(int k=1;k<=D;k++)
            {
                P[j].v[k]=w*P[j].v[k]+c1*URAND*(P[j].pbest[k]-P[j].p[k])+
                        c2*URAND*(gbest.p[k]-P[j].p[k]);
                if(P[j].v[k]>100)P[j].v[k]=100.0;
             //   printf("%f\n",P[j].v[k]);
                if(P[j].v[k]<-100)P[j].v[k]=-100.0;
                P[j].p[k]+=P[j].v[k];
                if(P[j].p[k]>100)P[j].p[k]=100;
            //    printf("%f\n",P[j].p[k]);
                if(P[j].p[k]<-100)P[j].p[k]=-100;
            }
            P[j].val=calc_val(P[j]);
            if(P[j].val<P[j].fitness)
            {
                for(int k=1;k<=D;k++)
                P[j].pbest[k]=P[j].p[k];

                P[j].fitness=P[j].val;
            }
            //收敛速度更快？？？？
            if(P[j].fitness<gbest.fitness)gbest=P[j];
        }
     //   w=0.7-times/G*0.5;
        //标准算法是在每代开始更新gbest的
        //for(int j=1;j<=P_num;j++)if(P[j].fitness<gbest.fitness)gbest=P[j];
    }
    gbest.val=calc_val(gbest);
    cout<<"the best value:  "<<gbest.val<<endl<<endl;
    printf("time:   %.6f\n\n\n\n",(double)clock()/CLOCKS_PER_SEC);
    system("pause");
    return 0;
}