粒子群算法源程序Word文档下载推荐.docx
《粒子群算法源程序Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《粒子群算法源程序Word文档下载推荐.docx（15页珍藏版）》请在冰豆网上搜索。
%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:
%[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）;
%开始进化，直到最大代数截至
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:
））;
%粒子速度更新（a）
%判断v的大小，限制v的绝对值小于
20
fork=1:
ifabs（v（j,k））>
20
rand（'
v（j,k）=20*rand（）;
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,:
%更新个体的历史极值
%以下更新全局的极值
%如果是最小的
[bestbest_index]=max（arr_present（:
值为min,相反则为max
ifbest>
gbest（end）&
（Region」n（arr_present（best_index,:
）,region））%女口果当前最好的结果
比以前的好，则更新最优值gbest,如果是最小的值为小于号，相反则为大于号
全局的极值
%混沌
xlhd=gbest（1:
if
（1）
forp=1:
25%次数
%1生成
cxl=rand（1,part_size）;
ifcxl（j）==0
cxl（j）=0.1;
ifcxl（j）==0.25
cxl（j）=0.26;
ifcxl（j）==0.5
cxl（j）=0.51;
ifcxl（j）==0.75
cxl（j）=0.76;
ifcxl（j）==1
cxl（j）=0.9;
%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）;
%4更新
cxl（j）=4*cxl（j）*（1-cxl（j））;
%当前代的最优粒子的适应度（取自）保存
best_record（:
i）=gbest;
%gbest:
一个行向量
pso=gbest;
%最优个体
display（gbest）;
figure;
plot（best_record（end,:
%最优解与代数的进化关系图
best=zeros（part_size,max_gen）;
part_size-1
best（i,:
）=best_record（i,:
pareto仁zeros（1,max_gen）;
pareto2=zeros（1,max_gen）;
pareto1（i）=f1（part_size,best（:
i））;
pareto2（i）=f2（part_size,best（:
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;
n-1
flag=flag&
（pos_present（1,j）>
=region（j,1））&
（pos_present（1,j）v=region（j,2））;
%初始化适应度
functionarr_fitness=ini_fit（pos_present,pop_size,part_size）
arr_fitness（k,1）=
fitness（part_size,pos_present（k,1:
%计算原始种群的适应
度
%%（***************************************************
***********************
%计算适应度
%%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;
「=叶（-10）*exp（（-0.2）*sqrt（x（iF2+x（i+1F2））;
%目标函数2
functionr=f2（n,x）
n
「=叶（abs（x（i））F0.8+5*sin（x（『3）;
%约束函数1
functionr=g1（n,x）
%约束函数2
functionr=g2（n,x）
%灰色关联度计算函数（越大相似性越好）
贷检（一组值）
应该取大于呢还是小于
%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））;
tn
if（abs（y0（i）-y1（i））>
m2）%
m2=abs（y0（i）-y1（i））;
kesa（i）=（m1+rou*m2）/（abs（y0（i）-y1（i））+rou*m2）;
gama=gama+kesa（i）;
gama=gama/tn;
%可行解的判决函数gn为约束条件的个数（暂时未用）n为解（粒子）的维数
functionbool=feasible（x,n）
%fori=1:
gn
r=max（0,g1（n,x）,g2（n,x））;
%判断约束条件
%end
if（r>
0）
bool=0;
%不可行解
else
bool=1;
%可行解
PSO粒子群算法解决旅行商问题的MATLAB源码
%粒子群算法求解旅行商问题
%ByIReij
closeall;
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）;
PopSize
Group（:
i）=randperm（CityNum）'
;
Group=Arrange（Group）;
%初始化粒子速度（即交换序）
Velocity=zeros（CityNum,PopSize）;
Velocity（:
i）=round（rand（1,CityNum）'
*CityNum）;
%round取整
%计算每个城市之间的距离
CityBetweenDistance=zeros（CityNum,CityNum）;
CityNum
CityBetweenDistance（i,j）=sqrt（（node（i,1）-node（j,1））A2+（node（i,2）-node（j,2））A2）;
end
%计算每条路径的距离
EachPathDis（i）=PathDistance（Group（:
i）'
CityBetweenDistanee）;
IndivdualBest=Group;
%记录各粒子的个体极值点位置，即个体找到的最短路径
IndivdualBestFitness=EachPathDis;
%记录最佳适应度值，即个体找到的最短路径的长度
[GlobalBestFitness,index]=min（EachPathDis）;
%找岀全局最优值和相应序号
%初始随机解
subplot（2,2,1）;
PathPlot（node,CityNum,index,IndivdualBest）;
随机解'
%寻优
while（lsStop==0）&
（Iteration<
Maxlteration）
%迭代次数递增
Iteration=Iteration+1;
%更新全局极值点位置，这里指路径
GlobalBest（:
i）=Group（:
index）;
%求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%更新各路径总距离
IsChange=EachPathDis<
IndivdualBestFitness;
%更新后的距离优于更新前的，记录序号
IndivdualBest（:
find（IsChange））=Group（:
find（IsChange））;
%更新个体最佳路径
更新个
IndivdualBestFitness=IndivdualBestFitness.*（~IsChange）+EachPathDis.*lsChange;
%体最佳路径距离
[GlobalBestFitness,index]=min（EachPathDis）;
%更新全局最佳路径，记录相应的序号
ifGlobalBestFitness==OldBestFitness%比较更新前和更新后的适应度值；
Num=Num+1;
%相等时记录加一；
OldBestFitness=GlobalBestFitness;
%不相等时更新适应度值，并记录清零
ifNum>
=20%多次迭代的适应度值相近时程序停止
lsStop=1;
BestFitness（lteration）=GlobalBestFitness;
%每一代的最优适应度
%最优解
subplot（2,2,2）;
优化解'
%进化曲线
subplot（2,2,3）;
plot（（1:
Iteration）,BestFitness（1:
Iteration））;
进化曲线'
%最小路径值
GlobalBestFitness
functionGroup=Arrange（Group）
[xy]=size（Group）;
[N01,index]=min（Group'
[],2）;
%找到最小值1
y
pop=Group（:
i）;
temp1=pop（[1:
index（i）_1]）;
temp2=pop（[index（i）:
x]）;
i）=[temp2'
tempi'
]'
endfunctionChangeNums=GenerateChangeNums（Group,BestVar）;
ChangeNums=zeros（x,y）;
pop=BestVar（:
%从BestVar取岀一个顺序
pop1=Group（:
%从粒子群中取岀对应的顺序
x%从BestVar的顺序中取岀一个序号
NoFromBestVar=pop（j）;
x%从对应的粒子顺序中取岀一个序号
NoFromGroup=pop1（k）;
if（NoFromBestVar==NoFromGroup）&
&
（j~=k）%两序号同且不在同一位置
ChangeNums（j,i）=k;
%交换子
pop1（k）=pop1（j）;
pop1（j）=NoFromGroup;
functionHold=HoldByOdds（Hold,Odds）
[x,y]=size（Hold）;
x
升级会员 