粒子群算法_matlab

最近在搞研究生数学建模，其中设涉及到一些很有意思的算法，贴出来与大家分享。

粒子群算法：

粒子群算法的思想就是模仿群鸟觅食的过程，每只鸟初始状态都是处于随机位置，且飞翔的方向也是随机的。他们都不知道食物在哪里，但是随着时间的推移，这些初始处于随机位置的鸟通过群内互相学习，信息共享和个体不断积累自身寻觅食物的经验，自组织集聚成一个群落，并逐渐朝唯一的目标——食物，前进。

在群鸟觅食模型中，每个个体可以被看成一个粒子，则鸟群可以被看成一个粒子群。其中第i个粒子的位置为Xi={x_i¹ ,x_i^2...x_i^d}，粒子的速度为Vi={v_i¹，v_i²...v_i^d}，d表示搜索空间的维度。通过每个位置计算出粒子群每个个体的适应值，根据适应值的大小确定其优劣势，并通过适应值来记录粒子个体经历过的最好位置P_i，记录整个粒子群经历过的最好位置P_g。粒子群采用下列算法对粒子群位置不断更新

以下为算法主要步骤

（1）初始化粒子群（速度和位置）、惯性因子、加速常数、最大迭代次数和算法终止的允许误差。

（2）评价每个粒子的初始适应值。

（3）将初始适应值作为每个粒子的局部最优适应值，将对应位置最为局部最优解。

（4）将最佳初始适应值作为全局最优适应值，将对应位置作为全局最优解。

（5）根据公式更新每个粒子的速度，并做限幅处理。

（6）根据公式更新粒子位置。

（7）比较当前适应值是否好于个体历史最优适应值，若好于则更新适应值和对应位置。

（8）在当前群中找到全局最优适应值，并将对应位置作为当前的粒子群全局最优解。

（9）重复（5）~（8），知道满足条件退出。

（10）输出全局最优适应值和其对应的位置。

function main()
%计算程序运行时间
tic;
%允许误差
e0=0.001; 
%最大迭代次数
maxnum=100;
%自变量个数
narvs=1;
%粒子群个体数
particlesize=30;
%加速度常数
c1=2;
c2=2;
%惯性因子
w=0.6;
%最大速度
vmax=0.8;
%初始化粒子，其在[-5,5]上随机分布
x=-5+10*rand(particlesize,narvs);
%初始化粒子的速度
v=2*rand(particlesize,narvs);
%适应值计算函数
fitness=inline('1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');
%计算初始的各个粒子的适应值
for i=1:particlesize
    for j=1:narvs
        f(i)=fitness(x(i,j));
    end
end
%初始的适应值最为个体最优适应值，对应的位置为个体最优位置
personalbest_x=x;
personalbest_faval=f;
%找出全局最好的适应值并获得其对应位置
[globalbest_faval i]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
%从1开始更新到最大迭代次数
k=1;
while k<=maxnum
    %根据新的位置获得适应值并更新个体历史最优适应值和对应位置
    for i=1:particlesize
        for j=i:narvs
            f(i)=fitness(x(i,j));
        end
        if f(i)<personalbest_faval(i)
            personalbest_faval(i)=f(i);
            personalbest_x(i,:)=x(i,:);
        end
    end
    [globalbest_faval i]=min(personalbest_faval);
    globalbest_x=personalbest_x(i,:);
    %根据公式计算出新的位置
    for i=1:particlesize
        v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))+c2*rand*(globalbest_x-x(i,:));
        %限制速度
        for j=1:narvs
            if v(i,j)>vmax;
                v(i,j)=vmax;
            elseif v(i,j)<-vmax;
                v(i,j)=-vmax;
            end
        end
        x(i,:)=x(i,:)+v(i,:);
    end
    if abs(globalbest_faval)<e0,break,end
    k=k+1;
end;
%绘图
value1=1/globalbest_faval-1;
value1=num2str(value1);
disp(strcat('the maximum value','=',value1));
value2=globalbest_x;
value2=num2str(value2);
disp(strcat('the corresponding coordinat','=',value2));
x=-5:0.01:5;
y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);
plot(x,y,'m-','linewidth',3);
hold on;
plot(globalbest_x,1/globalbest_faval-1,'kp','linewidth',4);
toc;

运行结果

猜你喜欢