%标准粒群优化算法程序
% 2007.1.9 By jxy
%测试函数:f(x,y)=100(x^2-y)^2+(1-x)^2,-2.048<x,y<2.048
%求解函数最小值
globalpopsize;%种群规模
%globalpopnum;%种群数量
globalpop;%种群
%globalc0;%速度惯性系数,为0—1的随机数
globalc1;%个体最优导向系数
globalc2;%全局最优导向系数
globalgbest_x;%全局最优解x轴坐标
globalgbest_y;%全局最优解y轴坐标
globalbest_fitness;%最优解
global best_in_history; %最优解变化轨迹
globalx_min;%x的下限
globalx_max;%x的上限
globaly_min;%y的下限
globaly_max;%y的上限
globalgen;%迭代次数
globalexetime;%当前迭代次数
globalmax_velocity;%最大速度
initial;%初始化
for exetime=1:gen
outputdata;%实时输出结果
adapting;%计算适应值
errorcompute(); %计算当前种群适值标准差
updatepop;%更新粒子位置
pause(0.01);
end
clear i;
clear exetime;
clear x_max;
clear x_min;
clear y_min;
clear y_max;
%程序初始化
gen=100;%设置进化代数
popsize=30;%设置种群规模大小
best_in_history(gen)=inf;%初始化全局历史最优解
best_in_history(=inf;%初始化全局历史最优解
max_velocity=0.3;%最大速度限制
best_fitness=inf;
%popnum=1;%设置种群数量
pop(popsize,8)=0;%初始化种群,创建popsize行5列的0矩阵
%种群数组第1列为x轴坐标,第2列为y轴坐标,第3列为x轴速度分量,第4列为y轴速度分量
%第5列为个体最优位置的x轴坐标,第6列为个体最优位置的y轴坐标
%第7列为个体最优适值,第8列为当前个体适应值
for i=1:popsize
pop(i,1)=4*rand()-2;%初始化种群中的粒子位置,值为-2—2,步长为其速度
pop(i,2)=4*rand()-2;%初始化种群中的粒子位置,值为-2—2,步长为其速度
pop(i,5)=pop(i,1); %初始状态下个体最优值等于初始位置
pop(i,6)=pop(i,2); %初始状态下个体最优值等于初始位置
pop(i,3)=rand()*0.02-0.01;%初始化种群微粒速度,值为-0.01—0.01,间隔为0.0001
pop(i,4)=rand()*0.02-0.01;%初始化种群微粒速度,值为-0.01—0.01,间隔为0.0001
pop(i,7)=inf;
pop(i,8)=inf;
end
c1=2;
c2=2;
x_min=-2;
y_min=- 2;
x_max=2;
y_max=2;
gbest_x=pop(1,1);%全局最优初始值为种群第一个粒子的位置
gbest_y=pop(1,2);
%适值计算
% 测试函数为f(x,y)=100(x^2-y)^2+(1-x)^2,-2.048<x,y<2.048
%计算适应值并赋值
for i=1:popsize
pop(i,8)=100*(pop(i,1)^2-pop(i,2))^2+(1-pop(i,1))^2;
ifpop(i,7)>pop(i,8)%若当前适应值优于个体最优值,则进行个体最优信息的更新
pop(i,7)=pop(i,8);%适值更新
pop(i,5:6)=pop(i,1:2);%位置坐标更新
end
end
%计算完适应值后寻找当前全局最优位置并记录其坐标
if best_fitness>min(pop(:,7))
best_fitness=min(pop(:,7));%全局最优值
gbest_x=pop(find(pop(:,7)==min(pop(:,7))),1);%全局最优粒子的位置
gbest_y=pop(find(pop(:,7)==min(pop(:,7))),2);
end
best_in_history(exetime)=best_fitness;%记录当前全局最优
%实时输出结果
%输出当前种群中粒子位置
subplot(1,2,1);
for i=1:popsize
plot(pop(i,1),pop(i,2),'b*');
holdon;
end
plot(gbest_x,gbest_y,'r.','markersize',20);axis([-2,2,-2,2]);
hold off;
subplot(1,2,2);
axis([0,gen,-0.00005,0.00005]);
if exetime-1>0
line([exetime-1,exetime],[best_in_history(exetime-1),best_fitness]);
hold on;
end
%粒子群速度与位置更新
%更新粒子速度
for i=1:popsize
pop(i,3)=rand()*pop(i,3)+c1*rand()*(pop(i,5)-pop(i,1))+c2*rand()*(gbest_x-pop(i,1));%更新速度
pop(i,4)=rand()*pop(i,4)+c1*rand()*(pop(i,6)-pop(i,2))+c2*rand()*(gbest_x-pop(i,2));
ifabs(pop(i,3))>max_velocity
if pop(i,3)>0
pop(i,3)=max_velocity;
else
pop(i,3)=-max_velocity;
end
end
ifabs(pop(i,4))>max_velocity
if pop(i,4)>0
pop(i,4)=max_velocity;
else
pop(i,4)=-max_velocity;
end
end
end
%更新粒子位置
for i=1:popsize
pop(i,1)=pop(i,1)+pop(i,3);
pop(i,2)=pop(i,2)+pop(i,4);
end
PSO算法-----VC++源码
///////////////////////////////////////////////////////////////
////
//求S=x1^2+x2^2+...+xn^2xi在(-50,100)之间的最小值//
////
///////////////////////////////////////////////////////////////
#include<iostream.h>
#include"stdlib.h"
#include"time.h"
#include"math.h"
#include<fstream.h>//处理文件
#include<windows.h>//处理系统时间
//随机数定义
#definerdint(i) (rand() % (int)(i))
#definerdft() (float)((double)rdint(16384)/(16383.0))
#definernd(a,b) (rdint((int)(b)-(int)(a)+1)+(int)(a))
//宏定义
#definePOPSIZE 40
#defineDIMENSION 10
#defineMAXITER 1000
//全局变量定义
floatw=0.9;
floatc1=1.8;
floatc2=1.8;
floatVMAX=100;
floatVMIN=-100;
floatXMIN=-50;
floatXMAX=100;
floatP[DIMENSION];
floatPBEST;
structindi
{
float number[DIMENSION];
float best[DIMENSION];
float bestfitness;
float fitness;
float speed[DIMENSION];
}individual[POPSIZE];
voidinitiate(void);
voidcalculation(int number);
voidglobalbest(int number);
voidlocalbest(int number);
//程序初始化定义
voidinitiate()
{
int i,j;
for(i=0;i<POPSIZE;i++)
{
for(j=0;j<DIMENSION;j++)
{
individual[i].number[j]=rdft()*(XMAX-XMIN)+XMIN;
}
}
for(i=0;i<POPSIZE;i++)
{
for(j=0;j<DIMENSION;j++)
{
individual[i].speed[j]=(VMAX-VMIN)*rdft()+VMIN;
}
}
for(i=0;i<POPSIZE;i++)
{
for(j=0;j<DIMENSION;j++)
{
individual[i].best[j]=individual[i].number[j];
}
}
for(i=0;i<POPSIZE;i++)
{
calculation(i);
}
for(i=0;i<POPSIZE;i++)
{
individual[i].bestfitness=individual[i].fitness;
}
globalbest(0);
}
//微粒历史最优位置修改程序
voidlocalbest(int number)
{
int i;
if(individual[number].bestfitness>individual[number].fitness)
{
for(i=0;i<DIMENSION;i++)
individual[number].best[i]=individual[number].number[i];
}
individual[number].bestfitness=individual[number].fitness;
}
//种群历史最优位置修改程序
voidglobalbest(int number)
{
int i,j,flag;
float s=0;
if (number==0)
{
s=individual[0].fitness;
flag=0;
for(i=1;i<POPSIZE;i++)
{
if(individual[i].fitness<s)
{
s=individual[i].fitness;
flag=i;
}
}
for(i=0;i<DIMENSION;i++)
{
P[i]=individual[flag].number[i];
}
PBEST=individual[flag].fitness;
}
else
{
for(i=0;i<POPSIZE;i++)
{
if(individual[i].bestfitness<PBEST)
{
for(j=0;j<DIMENSION;j++)
{
P[j]=individual[i].best[j];
}
PBEST=individual[i].bestfitness;
}
}
}
}
//适应值函数计算程序
voidcalculation(int num)
{
double s=0.0;
for(int i=0;i<DIMENSION;i++)
{
s+=pow(individual[num].number[i],2);
}
individual[num].fitness=s;
}
//更新速度和位置
voidUpdateSpeedandPosition()
{
int i,j;
i=0;j=0;
for(i=0;i<POPSIZE;i++)
{
for(j=0;j<DIMENSION;j++)
{individual[i].speed[j]=w*individual[i].speed[j]+c1*rdft()*(individual[i].best[j]-individual[i].number[j])+c2*rdft()*(P[j]-individual[i].number[j]);
if(individual[i].speed[j]>VMAX)
{
individual[i].speed[j]=VMAX;
}
if(individual[i].speed[j]<-VMAX)
{
individual[i].speed[j]=-VMAX;
}
individual[i].number[j]=individual[i].number[j]+individual[i].speed[j];
if(individual[i].number[j]<-XMAX)
{
individual[i].number[j]=-XMAX;
}
if(individual[i].number[j]>XMAX)
{
individual[i].number[j]=XMAX;
}
}
}
}
//主程序
voidmain()
{
int i,j,total;
i=0;
j=0;
total=0;
ofstream fout("C:\data.txt");//将数据输入到C:data.txt
SYSTEMTIMEmytime;
GetLocalTime(&mytime);
fout<<"程序开始执行时间是:"<<mytime.wYear<<"年"<<mytime.wMonth<<"月"
<<mytime.wDay<<"日"<<mytime.wHour<<"时"<<mytime.wMinute<<"分"
<<mytime.wSecond<<"秒"<<mytime.wMilliseconds<<"毫秒"<<"今天是星期"<<mytime.wDayOfWeek<<endl;
initiate();//初始化每个粒子
for (i=0;i<MAXITER;i++)
{
w=0.9-i*0.5/MAXITER;//动态修改权重
for (j=0;j<POPSIZE;j++) //计算每个粒子的适应值
{
calculation(j);
}
for (j=0;j<POPSIZE;j++) //计算更新粒子个体历史最优解
{
localbest(j);
}
for (j=0;j<POPSIZE;j++) //更新僵局最优解保存到PBEST和P[DIMENSION]中
{
globalbest(j);
}
UpdateSpeedandPosition();
total++;
cout<<"第"<<total<<"次迭代结果是:"<<"PBEST="<<PBEST<<""<<endl;
fout<<total<<""<<PBEST<<endl;
}
GetLocalTime(&mytime);
fout<<"程序结束执行时间是:"<<mytime.wYear<<"年"<<mytime.wMonth<<"月"
<<mytime.wDay<<"日"<<mytime.wHour<<"时"<<mytime.wMinute<<"分"
<<mytime.wSecond<<"秒"<<mytime.wMilliseconds<<"毫秒"<<"今天是星期"<<mytime.wDayOfWeek<<endl;
cout<<"最优解是:"<<PBEST<<endl;
int seq;
float pbest;
ifstream fin;
fin.open("C:\data.txt");
for (i=0;i<MAXITER;i++)
{fin>>seq>>pbest;
cout<<seq<<""<<pbest<<endl;
}
fin.close();
}