10.1 粒子群算法的MATLAB實現（2）

**10.1.3 **經典應用

粒子群算法經常與其他算法混合使用。混合策略就是將其他進化算法、傳統優化算法或其他技術應用到PSO中，用于提高粒子多樣性、增強粒子的全局探索能力，或者提高局部開發能力、增強收斂速度與精度。

常用的粒子群混合方法基于免疫的粒子群算法。該算法是在免疫算法的基礎上采用粒子群優化對抗體群體進行更新，可以解決免疫算法收斂速度慢的缺點。

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

① 確定學習因子c1和c 2 、粒子（抗體）群體個數M。

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

③ 產生免疫記憶粒子（抗體）：計算當前粒子（抗體）群體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)

% 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變化圖