10.1 粒子群算法的MATLAB實(shí)現(xiàn)（2）

**10.1.3 **經(jīng)典應(yīng)用

粒子群算法經(jīng)常與其他算法混合使用。混合策略就是將其他進(jìn)化算法、傳統(tǒng)優(yōu)化算法或其他技術(shù)應(yīng)用到PSO中，用于提高粒子多樣性、增強(qiáng)粒子的全局探索能力，或者提高局部開(kāi)發(fā)能力、增強(qiáng)收斂速度與精度。

常用的粒子群混合方法基于免疫的粒子群算法。該算法是在免疫算法的基礎(chǔ)上采用粒子群優(yōu)化對(duì)抗體群體進(jìn)行更新，可以解決免疫算法收斂速度慢的缺點(diǎn)。

基于免疫的混合粒子群算法步驟如下所示。

① 確定學(xué)習(xí)因子c1和c 2 、粒子（抗體）群體個(gè)數(shù)M。

② 由logistic回歸分析映射產(chǎn)生M個(gè)粒子（抗體）xi及其速度v i ，其中i=1,…,N，最后形成初始粒子（抗體）群體P 0 。

③ 產(chǎn)生免疫記憶粒子（抗體）：計(jì)算當(dāng)前粒子（抗體）群體P中粒子（抗體）的適應(yīng)度值并判斷算法是否滿足結(jié)束條件，如果滿足則結(jié)束并輸出結(jié)果，否則繼續(xù)運(yùn)行。

④ 更新局部和全局最優(yōu)解，并根據(jù)下面的公式更新粒子位置和速度。

⑤ 由logistic映射產(chǎn)生N個(gè)新的粒子（抗體）。

⑥ 基于濃度的粒子（抗體）選擇：用群體中相似抗體百分比計(jì)算產(chǎn)生N + M個(gè)新粒子（抗體）的概率，依照概率大小選擇N個(gè)粒子（抗體）形成粒子（抗體）群體P。然后轉(zhuǎn)入第③步。

算法流程圖如圖10-5所示。

圖10-5 免疫粒子群算法流程圖

將實(shí)現(xiàn)自適應(yīng)權(quán)重的優(yōu)化函數(shù)命名為PSO_immu.m，在MATLAB中編寫實(shí)現(xiàn)以上步驟的代碼：

function [x, y, Result] = PSO_immu(func, N, c1, c2, w, MaxDT, D, eps, DS, replaceP, minD, Psum)

format long;

%%%%%% 給定初始化條件 %%%%%%

% c1 = 2; % 學(xué)習(xí)因子1

% c2 = 2; % 學(xué)習(xí)因子2

% w = 0.8; % 慣性權(quán)重

% MaxDT = 100; % 最大迭代次數(shù)

% D = 2; % 搜索空間維數(shù)（未知數(shù)個(gè)數(shù))

% N = 100; % 初始化群體個(gè)體數(shù)目

% eps = 10^(-10); % 設(shè)置精度(在已知最小值時(shí)用)

% DS = 8; % 每隔DS次循環(huán)就檢查最優(yōu)個(gè)體是否變優(yōu)

% replaceP = 0. 5; % 粒子的概率大于replaceP將被免疫替換

% minD = 1e-10; % 粒子間的最小距離

% Psum = 0; % 個(gè)體最佳的和

range = 100;

count = 0;

%%%%%% 初始化種群的個(gè)體 %%%%%%

for i = 1 : N

for j = 1 : D

    x(i,j) = -range + 2 * range  * rand;    %隨機(jī)初始化位置

    v(i,j) = randn;         % 隨機(jī)初始化速度

end

end

%%%%%% 先計(jì)算各個(gè)粒子的適應(yīng)度,并初始化Pi和 Pg %%%%%%

for i = 1 : N

p(i) = feval(func,x(i, : ));

y(i, :) = x(i, :);

end

pg = x(1, :); % Pg為全局最優(yōu)

for i = 2 : N

if feval(func, x(i, :)) < feval(func, pg)

    pg = x(i, :);

end

end

%%%%%% 主循環(huán)，按照公式依次迭代，直到滿足精度要求 %%%%%%

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)

% break;

% end

%%%%%% 開(kāi)始進(jìn)行免疫 %%%%%%

if t > DS

    if (mod(t, DS) == 0) && ((Pbest(t - DS + 1) - Pbest(t)) < 1e-020)

        % 如果連續(xù)DS代數(shù)，群體中的最優(yōu)沒(méi)有明顯變優(yōu)，則進(jìn)行免疫

        % 在函數(shù)測(cè)試的過(guò)程中發(fā)現(xiàn)，經(jīng)過(guò)一定代數(shù)的更新，個(gè)體最優(yōu)不完全相等，但變化極小

        for i = 1 : N       % 先計(jì)算出個(gè)體最優(yōu)的和

            Psum = Psum + p(i);

        end

        for i = 1 : N       % 免疫程序

            for j = 1 : N   % 計(jì)算每個(gè)個(gè)體與個(gè)體i的距離

                distance(j) = abs(p(j)-p(i));

            end

            num = 0;

            for j = 1 : N   % 計(jì)算與第i個(gè)個(gè)體距離小于minD的個(gè)數(shù)

                if distance(j) < minD

                    num = num + 1;

                end

            end

            PF(i) = p(N - i + 1) / Psum;    % 計(jì)算適應(yīng)度概率

            PD(i) = num / N;        % 計(jì)算個(gè)體濃度

           

            a=rand;         % 隨機(jī)生成計(jì)算替換概率的因子

            PR(i) = a * PF(i) + (1 - a) * PD(i);    % 計(jì)算替換概率

        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

%%%%%% 最后給出計(jì)算結(jié)果 %%%%%%

x = pg(1, 1);

y = pg(1, 2);

Result = feval(func, pg);

%%%%%% 算法結(jié)束 %%%%%%

function probabolity(N, i)

PF = p(N - i) / Psum; % 適應(yīng)度概率

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; % 個(gè)體濃度

PR = a * PF + (1 - a) * PD; % 替換概率

例10-2：使用基于模擬退火的混合粒子群算法，求解函數(shù)的最小值。其中-10≤ xi ≤10，粒子數(shù)為50，學(xué)習(xí)因子均為2，退火常數(shù)取0.6，迭代步數(shù)為1000。

首先建立目標(biāo)函數(shù)代碼：

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)

運(yùn)行結(jié)果如下：

xm =

-0.735230545314342

fv =

1.241562924919382

得到目標(biāo)函數(shù)取最小值時(shí)的自變量xm變化圖，如圖10-6所示。

圖10-6 目標(biāo)函數(shù)取最小值時(shí)的自變量xm變化圖