最近在搞研究生数学建模，其中设涉及到一些很有意思的算法，贴出来与大家分享。
粒子群算法：
粒子群算法的思想就是模仿群鸟觅食的过程，每只鸟初始状态都是处于随机位置，且飞翔的方向也是随机的。他们都不知道食物在哪里，但是随着时间的推移，这些初始处于随机位置的鸟通过群内互相学习，信息共享和个体不断积累自身寻觅食物的经验，自组织集聚成一个群落，并逐渐朝唯一的目标——食物，前进。
在群鸟觅食模型中，每个个体可以被看成一个粒子，则鸟群可以被看成一个粒子群。其中第i个粒子的位置为Xi={xi1 ,xi2...xid}，粒子的速度为Vi={vi1，vi2...vid}，d表示搜索空间的维度。通过每个位置计算出粒子群每个个体的适应值，根据适应值的大小确定其优劣势，并通过适应值来记录粒子个体经历过的最好位置Pi，记录整个粒子群经历过的最好位置Pg。粒子群采用下列算法对粒子群位置不断更新
以下为算法主要步骤
（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;
运行结果
