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; }