基于模拟退火的粒子群优化算法（Matlab实现）

混合粒子群算法

混合粒子群算法是指借鉴其它一些智能优化算法的思想而形成的粒子群算法。除了粒子群算法外，还有遗传算法、模拟退火算法以及神经网络等智能算法，这些算法是目前应用比较广泛的智能算法，每种智能算法都有其特点，因此自然而然就有了结合各种智能算法的优点而形成的混合智能算法。下面介绍基于模拟退火算法的粒子群优化算法。
模拟退火算法在搜索过程中具有概率突跳的能力，能够有效避免搜索过程陷入局部极小解。模拟退火算法在退火过程中不但接受好的解，而且还以一定的概率接受差的解，同时这种概率受到温度参数的控制，其大小随着温度的下降而减小。

算法步骤

基于模拟退火的粒子群算法的步骤如下

随机初始化种群中各微粒的位置和速度；
评价每个微粒的适应度，将当前各微粒的位置和适应值存储在各微粒的pip_{i}pi中，将所有pbest中适应值最优个体的位置和适应值存储于pgp_{g}pg中；
确定初始温度；
根据下式确定当前温度下各pip_{i}pi的适配值

TF(pi)=e−[f(pi)−f(pg)]/t∑i=1Ne−[f(pi)−f(pg)]/tTF(p_{i})=\frac{e^{-[f(p_{i})-f(p_{g})]/t}}{\sum_{i=1}^{N}e^{-[f(p_{i})-f(p_{g})]/t}} TF(pi)=∑i=1Ne−[f(pi)−f(pg)]/te−[f(pi)−f(pg)]/t

采用轮盘赌策略从所有pip_{i}pi中确定全局最优的某个替代pg′p_{g}^{'}pg′，然后根据下式更新各微粒的速度和位置：

vi,j(t+1)=φ{vi,j(t)+c1r1[pi,j−xi,j(t)]+c2r2[pg,j′−xi,j(t)]}xi,j(t+1)=xi,j(t)+vi,j(t+1)v_{i,j}(t+1)=\varphi{\{v_{i,j}(t)+c_{1}r_{1}[p_{i,j}-x_{i,j}(t)]+c_{2}r_{2}[p_{g,j}^{'}-x_{i,j}(t)]\}} \\ x_{i,j}(t+1)=x_{i,j}(t)+v_{i,j}(t+1)\\ vi,j(t+1)=φ{vi,j(t)+c1r1[pi,j−xi,j(t)]+c2r2[pg,j′−xi,j(t)]}xi,j(t+1)=xi,j(t)+vi,j(t+1)
其中φ=2∣2−C−C2−4C∣,C=c1+c2;\varphi=\frac{2}{\vert 2-C-\sqrt{C^2-4C} \vert},C=c_1+c_2;φ=∣2−C−C2−4C∣2,C=c1+c2;

计算各微粒新的目标值，更新各微粒的pip_{i}pi值及群体的pgp_{g}pg值；
进行退温操作；
若满足停止条件（通常为预设的运算精度或迭代次数），搜索停止，输出结果，否则转step4；
初始温度和退温方式对算法有一定的影响，一般采用如下的初温和退温方式：
tk+1=λtk,t0=f(pg)/ln⁡5t_{k+1}=\lambda t_k,t_0=f(p_g)/\ln 5 tk+1=λtk,t0=f(pg)/ln5

算法的Matlab实现

function [xm,fv]=SimuAPSO(fitness,N,c1,c2,lamda,M,D)
% 功能：用基于模拟退火的粒子群优化算法求解无约束优化问题
% 待优化的目标函数：fitness
% 粒子数目：N
% 学习因子1：c1
% 学习因子2：c2
% 退火常数：lamda
% 最大迭代次数：M
% 自变量的个数：D
% 目标函数取最小值时的自变量值：xm
% 目标函数的最小值：fv
format long;
for i=1:Nfor j=1:Dx(i,j)=randn;   %随机初始化位置v(i,j)=randn;   %随机初始化速度end
end
for i=1:Np(i)=fitness(x(i,:));y(i,:)=x(i,:);
end
pg=x(N,:);      %pg为全局最优
for i=1:(N-1)if fitness(x(i,:))<fitness(pg)pg=x(i,:);end
end
T=fitness(pg)/log(5);       %初始温度
for t=1:MgroupFit=fitness(pg);for i=1:N               %当前温度下各个Pi的适应值Tfit(i)=exp(-(p(i)-groupFit)/T);endSumTfit=sum(Tfit);Tfit=Tfit/SumTfit;pBet=rand();for i=1:N               %用轮盘赌策略确定全局最优的某个替代值ComFit(i)=sum(Tfit(1:i));if pBet<=ComFit(i)pg_plus=x(i,:);break;endendC=c1+c2;ksi=2/abs(2-C-sqrt(C^2-4*C));       %速度压缩因子for i=1:Nv(i,:)=ksi*(v(i,:)+c1*randn*(y(i,:)-x(i,:))+c2*randn*(pg_plus-x(i,:)));x(i,:)=x(i,:)+v(i,:);if fitness(x(i,:))<p(i)p(i)=fitness(x(i,:));y(i,:)=x(i,:);endif p(i)<fitness(pg)pg=y(i,:);endendT=T*lamda;
end
xm=pg';
fv=fitness(pg);% [xm,fv]=SimuAPSO(@fitness,40,2.05,2.05,0.5,10000,5)

示例程序

求解如下最优化问题：
f(x)=[0.01+∑i=151i+(xi−1)2]−1,−10≤xi≤10,i=1,2,⋅⋅⋅,5f(x)=[0.01+\sum_{i=1}^{5} \frac{1}{i+(x_i-1)^2}]^{-1},-10 \leq x_i \leq 10,i=1,2,\sdot \sdot\sdot,5 f(x)=[0.01+i=1∑5i+(xi−1)21]−1,−10≤xi≤10,i=1,2,⋅⋅⋅,5
建立目标函数文件fitness.m

function F=fitness(x)
F=0;
for i=1:5F=F+1/(i+(x(i)-1)^2);
end
F=1/(0.01+F);
%以下是添加罚函数，将约束优化转变为无约束优化。
for i=1:5if x(i)>10||x(i)<-10F=inf;end
end

在Matlab命令行窗口输入：

 [xm,fv]=SimuAPSO(@fitness,40,2.05,2.05,0.5,10000,5)

所得结果为：

xm =1.0555919451626070.9639300049717460.6422443744351850.3994454832822640.270228145056420fv =0.447155413432555

参考文献

《精通MATLAB最优化计算》（第四版）龚纯、王正林，电子工业出版社