粒子群算法源程序.docx
《粒子群算法源程序.docx》由会员分享，可在线阅读，更多相关《粒子群算法源程序.docx（15页珍藏版）》请在冰豆网上搜索。
粒子群算法源程序
二维粒子群matlab源程序
%function[psoF]=pso_2D（）
%FUNCTIONPSOUSEParticleSwarmOptimization
Algorithm
%globalpresent;
%closeall;
clc;
clearall;
pop_size=10;%pop_size种群大小///粒子数量
part_size=2;%part_size粒子大小///粒子的维数
gbest=zeros（1,part_size+1）;%gbest当前搜索到的最小的值
max_gen=200;%max_gen最大迭代次数
%best=zeros（part_size,pop_size*part_size）;%xuan
region=zeros（part_size,2）;%设定搜索空间范围->解空间
region=10*[-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3];%
每一维设定不同范围（称之为解空间，不是可行域空间）
rand（'state',sum（100*clock））;%重置随机数发生器状态
%当前种群的信息矩阵，逐代进化的群体%当前位置,随机初始化
%一个10*3的随机的矩阵（初始化所有粒子的所有维数的位置值），其中最
后一列为
arr_present=ini_pos（pop_size,part_size）;
%初始化当前速度
%一个10*2的随机的矩阵（初始化所有粒子的所有维数的速度值）
v=ini_v（pop_size,part_size）;
%不是当前种群，可看作是一个外部的记忆体，存储每个粒子历史最优值（2维
数值）：
根据适应度更新！
%注意：
pbest数组10*3最后一列保存的是适应度
pbest=zeros（pop_size,part_size+1）;%pbest:
粒子以前搜索至U的最
优值，最后一列包括这些值的适应度
%1*80保存每代的最优值
best_record=zeros（part_size+1,max_gen）;%best_record
数组：
记录每一代的最好的粒子的适应度
w_max:
=0.9;
%w_max权系数最大值
w_min=
:
0.2;
%w_min权系数最小值
v_max=
:
2;
%最大速度,为粒子的范围宽度
c1=2;
%学习因子1
c2=2;
%学习因子2
%
%计算原始种群的适应度，及初始化
%
%注意：
传入的第一个参数是当前的粒子群体，ini_fit函数计算每个粒子的适
应度
%arr_present（:
end）是最后一列，保存每个粒子的适应值，是这样的！
xuan
arr_present（:
end）=ini_fit（arr_present,pop_size,part_size）;
%数组赋值，初始化每个粒子个体的历史最优值，以后会更新的
pbest=arr_present;%初始化各个粒子最优值
%找到当前群体中适应度最小的（在最后一列中寻找），best_value
%改为max，表示关联度最大
[best_valuebest_index]=max（arr_present（:
end））;%初始化全局最
优，即适应度为全局最小的值，根据需要也可以选取为最大值
%唯一的全局最优值，是当前代所有粒子中最好的一个
gbest=arr_present（best_index,:
）;
%因为是多目标，因此这个
%只是示意性的画出3维的
%x=[-3:
0.01:
3];
%y=[-3:
0.01:
3];
%[X,Y]=meshgrid（x,y）;
%Z1=（-10）*exp（（-0.2）*sqrt（XA2+YA2））;
%Z2=（abs（X））A0.8+abs（Y）A0.8+5*sin（XA3）+5*sin（丫八3）;
%z1=@（x,y）（-10）*exp（（-0.2）*sqrt（xA2+yA2））;
%z2=@（x,y）（abs（x））A0.8+abs（y）A0.8+5*sin（xA3）+5*sin（yA3）;%ezmeshc（z1）;gridon;
%ezmeshc（z2）;gridon;
%开始进化，直到最大代数截至
fori=1:
max_gen
%gridon;
%三维图象%多维图象是画不出来的
%ezmesh（z）,holdon,gridon;
%画出粒子群
%plot3（arr_present（:
1）,arr_present（:
2）,arr_present（:
3）,'*'）,h
oldoff;
%drawnow
%flush
%pause（0.01）;
线形递减权重
w=w_max-（w_max-w_min）*i/max_gen;%
%当前进化代数：
对于每个粒子进行更新和评价----->>>>>>>
forj=1:
pop_size
v（j,：
）=
w.*v（j,:
）+c1.*rand.*（pbest（j,1:
part_size）-arr_present（j,1:
part_size
））…
+c2.*rand*（gbest（1:
part_size）-arr_present（j,1:
part_size
））;%粒子速度更新（a）
%判断v的大小，限制v的绝对值小于
20
fork=1:
part_size
ifabs（v（j,k））>20
rand（'state',sum（100*clock））;
v（j,k）=20*rand（）;
end
end
%前几列是位置信息
arr_present（j,1:
part_size）=
arr_present（j,1:
part_size）+v（j,1:
part_size）;%粒子位置更新（b）
%最后一列是适应度
arr_present（j,end）=
fitness（part_size,arr_present（j,1:
part_size））;%适应度更新（保存至最
后一列）
%适应度评价与可行域限制
if
（arr_present（j,end）>pbest（j,end））&（Region」n（arr_present（j,:
）,regi
on））%根据条件更新pbest,如果是最小的值为小于号,相反则为大于号
pbest（j,:
）=arr_present（j,:
）;%更新个体的历史极值
end
end
%以下更新全局的极值
%如果是最小的
[bestbest_index]=max（arr_present（:
end））;值为min,相反则为max
ifbest>gbest（end）&
（Region」n（arr_present（best_index,:
）,region））%女口果当前最好的结果
比以前的好，则更新最优值gbest,如果是最小的值为小于号，相反则为大于号
全局的极值
gbest=arr_present（best_index,:
）;%
end
%混沌
xlhd=gbest（1:
part_size）;
if
（1）
forp=1:
25%次数
%1生成
cxl=rand（1,part_size）;
forj=1:
part_size
ifcxl（j）==0
cxl（j）=0.1;
end
ifcxl（j）==0.25
cxl（j）=0.26;
end
ifcxl（j）==0.5
cxl（j）=0.51;
end
ifcxl（j）==0.75
cxl（j）=0.76;
end
ifcxl（j）==1
cxl（j）=0.9;
end
end
%2映射
al=-30;bl=30;
rxl=al+（bl-al）*cxl;
%3搜索
bate=0.1;
xlhd=xlhd+bate*rxl;
iffitness（part_size,xlhd）>gbest（end）
gbest（1:
part_size）=xlhd;
gbest（end）=fitness（part_size,xlhd）;
end
%4更新
forj=1:
part_size
cxl（j）=4*cxl（j）*（1-cxl（j））;
end
end
end
%混沌
%当前代的最优粒子的适应度（取自）保存
best_record（:
i）=gbest;%gbest:
一个行向量
end
pso=gbest;%最优个体
display（gbest）;
figure;
plot（best_record（end,:
））;%最优解与代数的进化关系图
best=zeros（part_size,max_gen）;
fori=1:
part_size-1
best（i,:
）=best_record（i,:
）;
end
pareto仁zeros（1,max_gen）;
pareto2=zeros（1,max_gen）;
fori=1:
max_gen
pareto1（i）=f1（part_size,best（:
i））;
pareto2（i）=f2（part_size,best（:
i））;
end
figure;
i=1:
max_gen;
%plot（i,pareto1（i）,'r*',i,pareto2（i）,'g*'）;
plot（pareto1（i）,pareto2（i）,'r+'）;
xlabel（'f1'）;ylabel（'f2'）;
title（'Pareto曲线'）;
%figure;
%plot（,f2（best_record）,）;
%movie2avi（F,'pso_2D1.avi','compression','MSVC'）;
%子函数
%
%
%返回随机的位置
functionini_present=ini_pos（pop_size,part_size）
ini_present=10*3*rand（pop_size,part_size+1）;%初始化当前
粒子位置,使其随机的分布在工作空间
%返回一个随机的矩阵，10*（2+1），最后一列将用来保存适应度
%返回随机的速度
functionini_velocity=ini_v（pop_size,part_size）
ini_velocity=20*（rand（pop_size,part_size））;%初始化当前粒子速
度,使其随机的分布在速度范围内
%判断是否处于范围内
functionflag=Region_in（pos_present,region）
[mn]=size（pos_present）;%1*11n返回解的维数10
flag=1;
forj=1:
n-1
flag=flag&（pos_present（1,j）>=region（j,1））&
（pos_present（1,j）v=region（j,2））;
end
%初始化适应度
functionarr_fitness=ini_fit（pos_present,pop_size,part_size）
fork=1:
pop_size
arr_fitness（k,1）=
fitness（part_size,pos_present（k,1:
part_size））;%计算原始种群的适应
度
end
%%（***************************************************
***********************
%计算适应度
%%i***************************************************
***********************functionfit=fitness（n,xp）
%需要求极值的函数,本例即peaks函数
%y0=[-85.4974,-29.9217];%注意：
这是基准序列，也就是单个最优的极
值
y0=[-9.9907,-7.7507];
%y0=[-39.6162,-18.4561];
%y0=[-86.8312,-29.9217];
y1=[f1（n,xp）,f2（n,xp）];%n为粒子维数
fit=graydegree（2,y0,y1）;%关联度在某种意义上就是适应度
%目标函数1
functionr=f1（n,x）
r=0;
fori=1:
n-1
「=叶（-10）*exp（（-0.2）*sqrt（x（iF2+x（i+1F2））;
end
%目标函数2
functionr=f2（n,x）
r=0;
fori=1:
n
「=叶（abs（x（i））F0.8+5*sin（x（『3）;
end
%约束函数1
functionr=g1（n,x）
r=0;
fori=1:
n
r=0;
end
%约束函数2
functionr=g2（n,x）
r=0;
fori=1:
n
r=0;
end
%灰色关联度计算函数（越大相似性越好）
贷检（一组值）
应该取大于呢还是小于
%tn目标函数个数x0基准序列（一组值）x1
functiongama=graydegree（tn,y0,y1）
gama=0;
rou=0.5;
kesa=zeros（tn,1）;
m1=abs（y0
（1）-y1
（1））;
m2=abs（y0
（1）-y1
（1））;
fori=1:
tn
if（abs（y0（i）-y1（i））>m2）%
m2=abs（y0（i）-y1（i））;
end
end
fori=1:
tn
kesa（i）=（m1+rou*m2）/（abs（y0（i）-y1（i））+rou*m2）;
gama=gama+kesa（i）;
end
gama=gama/tn;
%可行解的判决函数gn为约束条件的个数（暂时未用）n为解（粒子）的维数
functionbool=feasible（x,n）
r=0;
%fori=1:
gn
r=max（0,g1（n,x）,g2（n,x））;%判断约束条件
%end
if（r>0）
bool=0;%不可行解
else
bool=1;%可行解
end
PSO粒子群算法解决旅行商问题的MATLAB源码
PSO粒子群算法解决旅行商问题的MATLAB源码
%粒子群算法求解旅行商问题
%ByIReij
closeall;
clearall;
PopSize=500;%种群大小
CityNum=14;%城市数
OldBestFitness=0;%旧的最优适应度值
Iteration=0;%迭代次数
Maxlteration=2000;%最大迭代次数
lsStop=0;%程序停止标志
Num=0;%取得相同适应度值的迭代次数
c1=0.5;%认知系数
c2=0.7;%社会学习系数
w=0.96-lteration/Maxlteration;%惯性系数，随迭代次数增加而递减
%节点坐标
node=[16.4796.10;16.4794.44;20.0992.54;22.3993.37;25.2397.24;...
22.0096.05;20.4797.02;17.2096.29;16.3097.38;14.0598.12;...
16.5397.38;21.5295.59;19.4197.13;20.0994.55];
%初始化各粒子，即产生路径种群
Group=ones（CityNum,PopSize）;
fori=1:
PopSize
Group（:
i）=randperm（CityNum）';
end
Group=Arrange（Group）;
%初始化粒子速度（即交换序）
Velocity=zeros（CityNum,PopSize）;
fori=1:
PopSize
Velocity（:
i）=round（rand（1,CityNum）'*CityNum）;%round取整
end
%计算每个城市之间的距离
CityBetweenDistance=zeros（CityNum,CityNum）;
fori=1:
CityNum
forj=1:
CityNum
CityBetweenDistance（i,j）=sqrt（（node（i,1）-node（j,1））A2+（node（i,2）-node（j,2））A2）;end
end
%计算每条路径的距离
fori=1:
PopSize
EachPathDis（i）=PathDistance（Group（:
i）',CityBetweenDistanee）;
end
IndivdualBest=Group;%记录各粒子的个体极值点位置，即个体找到的最短路径
IndivdualBestFitness=EachPathDis;%记录最佳适应度值，即个体找到的最短路径的长度
[GlobalBestFitness,index]=min（EachPathDis）;%找岀全局最优值和相应序号
%初始随机解
figure;
subplot（2,2,1）;
PathPlot（node,CityNum,index,IndivdualBest）;
title（'随机解'）;
%寻优
while（lsStop==0）&（Iteration%迭代次数递增
Iteration=Iteration+1;
%更新全局极值点位置，这里指路径
fori=1:
PopSize
GlobalBest（:
i）=Group（:
index）;
end
%求pij-xij,pgj-xij交换序，并以概率c1,c2的保留交换序
pij_xij=GenerateChangeNums（Group,IndivdualBest）;
pij_xij=HoldByOdds（pij_xij,c1）;
pgj_xij=GenerateChangeNums（Group,GlobalBest）;
Pgj_xij=HoldByOdds（pgj_xij,c2）;
%以概率w保留上一代交换序
Velocity=HoldByOdds（Velocity,w）;
Group=PathExchange（Group,Velocity）;%根据交换序进行路径交换
Group=PathExchange（Group,pij_xij）;
Group=PathExchange（Group,pgj_xij）;
fori=1:
PopSize%更新各路径总距离
EachPathDis（i）=PathDistance（Group（:
i）',CityBetweenDistanee）;
end
IsChange=EachPathDisIndivdualBest（:
find（IsChange））=Group（:
find（IsChange））;%更新个体最佳路径
更新个
IndivdualBestFitness=IndivdualBestFitness.*（~IsChange）+EachPathDis.*lsChange;%体最佳路径距离
[GlobalBestFitness,index]=min（EachPathDis）;%更新全局最佳路径，记录相应的序号
ifGlobalBestFitness==OldBestFitness%比较更新前和更新后的适应度值；
Num=Num+1;%相等时记录加一；
else
OldBestFitness=GlobalBestFitness;%不相等时更新适应度值，并记录清零
Num=0;
end
ifNum>=20%多次迭代的适应度值相近时程序停止
lsStop=1;
end
BestFitness（lteration）=GlobalBestFitness;%每一代的最优适应度
end
%最优解
subplot（2,2,2）;
PathPlot（node,CityNum,index,IndivdualBest）;
title（'优化解'）;
%进化曲线
subplot（2,2,3）;
plot（（1:
Iteration）,BestFitness（1:
Iteration））;
gridon;
title（'进化曲线'）;
%最小路径值
GlobalBestFitness
functionGroup=Arrange（Group）
[xy]=size（Group）;
[N01,index]=min（Group',[],2）;%找到最小值1
fori=1:
y
pop=Group（:
i）;
temp1=pop（[1:
index（i）_1]）;
temp2=pop（[index（i）:
x]）;
Group（:
i）=[temp2'tempi']';
endfunctionChangeNums=GenerateChangeNums（Group,BestVar）;
[xy]=size（Group）;
ChangeNums=zeros（x,y）;
fori=1:
y
pop=BestVar（:
i）;%从BestVar取岀一个顺序
pop1=Group（:
i）;%从粒子群中取岀对应的顺序
forj=1:
x%从BestVar的顺序中取岀一个序号
NoFromBestVar=pop（j）;
fork=1:
x%从对应的粒子顺序中取岀一个序号
NoFromGroup=pop1（k）;
if（NoFromBestVar==NoFromGroup）&&（j~=k）%两序号同且不在同一位置
ChangeNums（j,i）=k;%交换子
pop1（k）=pop1（j）;
pop1（j）=NoFromGroup;
end
end
end
end
functionHold=HoldByOdds（Hold,Odds）
[x,y]=size（Hold）;
fori=1:
x
升级会员 