单目标----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; %学习因⼦1c2=1.4962; %学习因⼦2w=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; %随机初始化速度 endend

%------先计算各个粒⼦的适应度，并初始化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)%------进⼊主要循环，按照公式依次迭代，直到满⾜精度要求------------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)if p(i)Pbest(t)=f(pg,D);end

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

disp('*************************************************************')disp('函数的全局最优位置为：')Solution=pg'

disp('最后得到的优化极值为：')Result=f(pg,D)

disp('*************************************************************')%------算法结束---DreamSun GL & HF-----------------------------------

最后再贴下我⾃⼰写的代码吧（测试函数依旧⽤的JADE⾥的第⼀个，精度也⾮常精确了~~）代码风格依旧很挫。。。。

#include#include#include#include//#include#include#includeusing 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\

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\

if(P[j].p[k]<-100)P[j].p[k]=-100; }

P[j].val=calc_val(P[j]); if(P[j].valfor(int k=1;k<=D;k++)

P[j].pbest[k]=P[j].p[k];

P[j].fitness=P[j].val; }

//收敛速度更快

if(P[j].fitness// w=0.7-times/G*0.5;

//标准算法是在每代开始更新gbest的

//for(int j=1;j<=P_num;j++)if(P[j].fitnessgbest.val=calc_val(gbest);

cout<<\"the best value: \"<printf(\"time: %.6f\\n\\n\\n\\n\",(double)clock()/CLOCKS_PER_SEC); system(\"pause\"); return 0;}