学习,是一个日积月累的过程……
 

粒子群优化算法(PSO)MATLAB源程序

作者:Aaron    博客:微亮的烛光

题记:由于最近有网友根据我之前编写的粒子群算法向我提问,就重新看了一下自己之前的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


 
评论(5)
 
热度(7)
© 微亮的烛光|Powered by LOFTER