粒子群算法源程序.docx
- 文档编号:7706910
- 上传时间:2023-01-25
- 格式:DOCX
- 页数:15
- 大小:20.87KB
粒子群算法源程序.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=EachPathDis IndivdualBest(: 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 粒子 算法 源程序