10.1 粒子群算法的MATLAB实现(2)
10.1.3 经典应用
粒子群算法经常与其他算法混合使用。混合策略就是将其他进化算法、传统优化算法或其他技术应用到PSO中,用于提高粒子多样性、增强粒子的全局探索能力,或者提高局部开发能力、增强收敛速度与精度。
常用的粒子群混合方法基于免疫的粒子群算法。该算法是在免疫算法的基础上采用粒子群优化对抗体群体进行更新,可以解决免疫算法收敛速度慢的缺点。
基于免疫的混合粒子群算法步骤如下所示。
①确定学习因子c1和c2、粒子(抗体)群体个数M。
②由logistic回归分析映射产生M个粒子(抗体)xi及其速度vi,其中i=1,…,N,最后形成初始粒子(抗体)群体P0。
③产生免疫记忆粒子(抗体):计算当前粒子(抗体)群体P中粒子(抗体)的适应度值并判断算法是否满足结束条件,如果满足则结束并输出结果,否则继续运行。
④更新局部和全局最优解,并根据下面的公式更新粒子位置和速度。
⑤由logistic映射产生N个新的粒子(抗体)。
⑥基于浓度的粒子(抗体)选择:用群体中相似抗体百分比计算产生N + M个新粒子(抗体)的概率,依照概率大小选择N个粒子(抗体)形成粒子(抗体)群体P。然后转入第③步。
算法流程图如图10-5所示。
图10-5 免疫粒子群算法流程图
将实现自适应权重的优化函数命名为PSO_immu.m,在MATLAB中编写实现以上步骤的代码:
function [x, y, Result] = PSO_immu(func, N, c1, c2, w, MaxDT, D, eps, DS, replaceP, minD, Psum) format long; %%%%%% 给定初始化条件 %%%%%% % c1 = 2; % 学习因子1 % c2 = 2; % 学习因子2 % w = 0.8; % 惯性权重 % MaxDT = 100; % 最大迭代次数 % D = 2; % 搜索空间维数(未知数个数) % N = 100; % 初始化群体个体数目 % eps = 10^(-10); % 设置精度(在已知最小值时用) % DS = 8; % 每隔DS次循环就检查最优个体是否变优 % replaceP = 0. 5; % 粒子的概率大于replaceP将被免疫替换 % minD = 1e-10; % 粒子间的最小距离 % Psum = 0; % 个体最佳的和 range = 100; count = 0; %%%%%% 初始化种群的个体 %%%%%% for i = 1 : N for j = 1 : D x(i,j) = -range + 2 * range * rand; %随机初始化位置 v(i,j) = randn; % 随机初始化速度 end end %%%%%% 先计算各个粒子的适应度,并初始化Pi和 Pg %%%%%% for i = 1 : N p(i) = feval(func,x(i, : )); y(i, :) = x(i, :); end pg = x(1, :); % Pg为全局最优 for i = 2 : N if feval(func, x(i, :)) < feval(func, pg) pg = x(i, :); end end %%%%%% 主循环,按照公式依次迭代,直到满足精度要求 %%%%%% for t = 1 : MaxDT for i = 1 : N v(i, :) = w * v(i, :) + c1 * rand * (y(i, :) - x(i, :) ) + c2 * rand * (pg - x (i, :)); x(i, :) = x(i, :) + v(i, :); if feval(func,x(i, :)) < p(i) p(i) = feval(func, x(i, :)); y(i, :) = x(i, :); end if p(i) < feval(func, pg) pg = y(i, :); subplot(1, 2, 1); bar(pg, 0.25); axis( [0 3 -40 40]); title ( [ 'Iteration', num2str(t)]); pause(0.1); subplot(1, 2, 2); plot(pg(1, 1), pg(1, 2), 'rs', 'MarkerFacecolor', 'r', 'Markersize', 8) hold on; plot(x(:, 1),x(:, 2), 'k.'); set(gca, 'Color', 'g') hold off; grid on; axis([-100 100 -100 100]); title(['Global Min = ' , num2str(p(i))]); xlabel(['Min_x= ', num2str(pg(1, 1)), 'Min_y= ', num2str(pg(1, 2))]); end end Pbest(t) = feval(func, pg); % if Foxhole(pg,D)<eps % 如果结果满足精度要求则跳出循环 % break; % end %%%%%% 开始进行免疫 %%%%%% if t > DS if (mod(t, DS) == 0) && ((Pbest(t - DS + 1) - Pbest(t)) < 1e-020) % 如果连续DS代数,群体中的最优没有明显变优,则进行免疫 % 在函数测试的过程中发现,经过一定代数的更新,个体最优不完全相等,但变化极小 for i = 1 : N % 先计算出个体最优的和 Psum = Psum + p(i); end for i = 1 : N % 免疫程序 for j = 1 : N % 计算每个个体与个体i的距离 distance(j) = abs(p(j)-p(i)); end num = 0; for j = 1 : N % 计算与第i个个体距离小于minD的个数 if distance(j) < minD num = num + 1; end end PF(i) = p(N - i + 1) / Psum; % 计算适应度概率 PD(i) = num / N; % 计算个体浓度 a=rand; % 随机生成计算替换概率的因子 PR(i) = a * PF(i) + (1 - a) * PD(i); % 计算替换概率 end for i = 1 : N if PR(i) > replaceP x(i, :) = -range + 2 * range * rand(1, D); count=count+1; end end end end end %%%%%% 最后给出计算结果 %%%%%% x = pg(1, 1); y = pg(1, 2); Result = feval(func, pg); %%%%%% 算法结束 %%%%%% function probabolity(N, i) PF = p(N - i) / Psum; % 适应度概率 disp(PF); for jj = 1 : N distance(jj) = abs(P(jj)-P(i)); end num = 0; for ii = 1 : N if distance(ii) < minD num = num + 1; end end PD = num/N; % 个体浓度 PR = a * PF + (1 - a) * PD; % 替换概率
例10-2:使用基于模拟退火的混合粒子群算法,求解函数的最小值。其中-10≤ xi≤10,粒子数为50,学习因子均为2,退火常数取0.6,迭代步数为1000。
首先建立目标函数代码:
function y = immuFunc(x) y = (cos(x(1)^2 + x(2)^2) - 1)/((1 + (x(1)^2 - x(2)^2))^2) + 0.5; end
在MATLAB命令行窗口中输入代码:
[xm, fv] = PSO_immu(@immuFunc, 50, 2, 2, 0.8, 100, 5, 0.0000001, 10, 0.6, 0.0000000000000000001, 0)
运行结果如下:
xm = -0.735230545314342 fv = 1.241562924919382
得到目标函数取最小值时的自变量xm变化图,如图10-6所示。
图10-6 目标函数取最小值时的自变量xm变化图