粒子群算法源程序Word文档下载推荐.docx

粒子群算法源程序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