粒子群优化算法(PSO)的MATLAB源程序
题记:由于最近有网友根据我之前编写的粒子群算法向我提问,就重新看了一下自己之前的MATLAB程序,发现自己当时编写的程序从易读性和简洁性上都惨不忍睹,遂依据自己对MATLAB编程最新的理解,以矩阵思想重新编写了PSO程序,极大的简化了程序。同时在征得网友同意后将一个算例分享出来,希望可以帮助更多的人。
说明:粒子群优化算法(Particle Swarm Optimization, PSO)是群智能优化算法之一,具有便于实现和收敛速度快等优点。本人在研究这个算法的时候,编写了一些测试的MATLAB源程序,在此分享,以供学习交流之用。这个是最基本的粒子群优化算法。
文件结构:
主函数:pso.m
适应度函数:fitness.m
源代码:
下载地址:pso.rar
---------------------------pso.m---------------------------
function [ ] = pso( )
%PSO Algorithm
%求解变量需要归一化在(0,1)范围
%求解目标函数为使适应度最小
clear all;
%设置数据格式
format short;
%设置粒子群算法参数
p_num=50; %粒子数量(种群规模)
p_size=2; %粒子维数(求解变量个数)
i_max=50; %最大迭代次数
%%%%%%%%初始化迭代结果存储矩阵
gbest=zeros(p_size,i_max); %全局最优粒子位置矩阵
gbestf=zeros(1,i_max); %全局最优粒子适应度矩阵
%%%%%%%%初始化粒子群
i=1; %迭代次数,index ofiteration
particle=rand(p_size,p_num); %粒子位置(0,1)
particle=(ones(p_size,1)*(sum(particle)<1)).*particle;%目标函数额外限制条件p1+p2<1,否则至0
pfitness=fitness(particle); %粒子适应度
pbest=particle; %局部最优粒子
pbestf=pfitness; %局部最优粒子适应度
[gbestf(1,i),gIndex]=min(pbestf); %全局最优粒子适应度及索引
gbest(:,i)=pbest(:,gIndex); %全局最优粒子位置
velocity=-1+2.*rand(p_size,p_num); %粒子速度(-1,1)
%%%%%%%%迭代计算
for i=2:i_max
%粒子位置 - 位置更新公式&取值范围限制条件(0,1)
p=particle+velocity;
p=p.*(p>=0).*(p<=1)+(p<0).*(1e-4)+(p>1).*(9e-4);
p=(ones(p_size,1)*(sum(p)<1)).*p+(ones(p_size,1)*(sum(p)>=1)).*particle; %目标函数额外限制条件p1+p2<1
%粒子适应度 - 目标函数
pf=fitness(p);
%局部最优粒子 - 取适应度小的粒子
pbest=(ones(p_size,1)*(pf<=pfitness)).*p+(ones(p_size,1)*(pf>pfitness)).*particle;
%局部最优粒子适应度
pbestf=(pf<=pfitness).*pf+(pf>pfitness).*pfitness;
%全局最优粒子适应度及索引 - 取适应度全局最小的粒子
[gbestf(1,i),gIndex]=min(pbestf);
%全局最优粒子位置
gbest(:,i)=pbest(:,gIndex);
%粒子速度 - 速度更新公式&取值范围限制条件(-1,1)
v=velocity+1.3*rand*(pbest-p)+1.3*rand*(gbest(:,i)*ones(1,p_num)-p);
v=v.*(v>=-1).*(v<=1)+(v<-1).*(-9e-4)+(v>1).*(9e-4);
%更新粒子位置
particle=pbest;
%更新粒子速度
velocity=v;
end
save('ResultGbest','gbest','gbestf','-ascii','-tabs','-append');
figure=plot(gbestf);
xlabel('迭代次数');
title('全局最优粒子适应度');
saveas(figure,'ResultGbest.fig');
saveas(figure,'ResultGbest.bmp');
end
--------------------------fitness.m--------------------------
function [ f ] =fitness( p )
%The TragetFunction for PSO to calculate the fitness
[p_size,p_num]=size(p);
f=100*(30*ones(1,p_num)-10*p(1,:)-24*p(2,:));
end
算例:
数学规划模型min(x)=10*x1+3*x2+15*x3,
约束条件是x1+x2+x3=200,x1,x2,x3都大于零
可以理解为要求得使min(x)最小的x1, x2, x3的取值
求解:
根据x1+x2+x3=200,可以将问题降维
取x3=200-x1-x2,且x3>0,
有min(x)=10*x1+3*x2+15*(200-x1-x2)=3000-5*x1-12*x2,
约束条件为200>x1>0,200>x2>0,x1+x2<200,
归一化使p1=x1/200, p2=x2/200, 有p1+p2<1,
min(p)=3000-1000*p1-2400*p2=100*(30-10*p1-24*p2),
即只要求得使min(p)最小的p1和p2即可
说明:
在原程序中的变量对应关系
i:迭代次数
particle(1,:):即求解问题中的p1
particle(2,:):即求解问题中的p2
fitness():即求解的目标函数min(p)
所求的结果存于gbest矩阵中,每迭代一次存储一次结果
程序的所有运行结果存于文件ResultGbest.mat
结果:
PSO的近似结果为:
min(p)=602
p1=0.999,p2=0.0001即x1=199.8, x2=0.02, x3=0.18
寻优结果:
本文是原创文章,转载请保留原作者和出处信息。
由于个人水平有限,不足之处还望多多包涵,欢迎批评指正。
By Aaron