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

粒子群算法（Particle Swarm Optimization，PSO）屬于進(jìn)化算法的一種，和模擬退火算法相似，它也是從隨機(jī)解出發(fā)，通過迭代尋找最優(yōu)解。它也是通過適應(yīng)度來評價解的品質(zhì)，但它比遺傳算法規(guī)則更為簡單，沒有遺傳算法的“交叉”（Crossover）和“變異”（Mutation）操作，它通過追隨當(dāng)前搜索到的最優(yōu)值來尋找全局最優(yōu)。

**10.1.1 **基本原理

PSO可以用于解決優(yōu)化問題。在PSO中，每個優(yōu)化問題的潛在解都是搜索空間中的一只鳥，稱為粒子。所有的粒子都有一個由被優(yōu)化的函數(shù)決定的適值（Fitness Value），每個粒子還有一個速度決定它們“飛行”的方向和距離。然后粒子們就追隨當(dāng)前的最優(yōu)粒子在解空間中搜索。

粒子位置的更新方式如圖10-1所示。

圖10-1 粒子位置的更新方式

其中，x表示粒子起始位置，v表示粒子“飛行”的速度，p表示搜索到的粒子的最優(yōu)位置。

PSO初始化為一群隨機(jī)粒子（隨機(jī)解），然后通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個極值來更新自己：一個是粒子本身所找到的最優(yōu)解，這個解稱為個體極值；另一個極值是整個種群目前找到的最優(yōu)解，這個極值是全局極值。

另外，也可以不用整個種群而只用其中一部分作為粒子的鄰居，那么在所有鄰居中的極值就是局部極值。

假設(shè)在一個D維的目標(biāo)搜索空間中，有N個粒子組成一個群落，其中第i個粒子表示為一個D維的向量

第i個粒子的“飛行”速度也是一個D維的向量，記為

第i個粒子迄今為止搜索到的最優(yōu)位置稱為個體極值，記為

整個粒子群迄今為止搜索到的最優(yōu)位置為全局極值，記為

在找到這兩個最優(yōu)值時，粒子根據(jù)如下公式來更新自己的速度和位置：

其中，c1和c2為學(xué)習(xí)因子，也稱加速常數(shù)（Acceleration Constant）；r1和r2為[0,1]范圍內(nèi)的均勻隨機(jī)數(shù)。

式右邊由三部分組成：

● 第一部分為“慣性”（Inertia）或“動量”（Momentum）部分，反映了粒子的運(yùn)動“習(xí)慣（Habit）”，代表粒子有維持自己先前速度的趨勢。

● 第二部分為“認(rèn)知”（Cognition）部分，反映了粒子對自身歷史經(jīng)驗(yàn)的記憶或回憶，代表粒子有向自身歷史最佳位置逼近的趨勢。

● 第三部分為“社會”（Social）部分，反映了粒子間協(xié)同合作與知識共享的群體歷史經(jīng)驗(yàn)，代表粒子有向群體或鄰域歷史最佳位置逼近的趨勢。

由于粒子群算法具有高效的搜索能力，因此有利于得到多目標(biāo)意義下的最優(yōu)解；通過代表整個解集種群，按并行方式同時搜索多個非劣解，即搜索到多個Pareto最優(yōu)解。

同時，粒子群算法的通用性比較好，適合處理多種類型的目標(biāo)函數(shù)和約束，并且容易與傳統(tǒng)的優(yōu)化方法結(jié)合，從而改進(jìn)自身的局限性，更高效地解決問題。因此，將粒子群算法應(yīng)用于解決多目標(biāo)優(yōu)化問題上具有很大的優(yōu)勢。

**10.1.2 **程序設(shè)計

基本粒子群算法的流程圖如圖10-2所示。其具體過程如下：

圖10-2 基本粒子群算法流程圖

① 初始化粒子群，包括群體規(guī)模N、每個粒子的位置xi和速度v i 。

② 計算每個粒子的適應(yīng)度值F it [i]。

③ 對每個粒子，用它的適應(yīng)度值F it [i]和個體極值p best (i)比較，如果F it [i] > p best (i)，則用F it [i]替換p best (i)。

④ 對每個粒子，用它的適應(yīng)度值F it [i]和個體極值g best (i)比較，如果F it [i] > p best (i)，則用F it [i]替換g best (i)。

⑤ 更新粒子的速度vi和位置x i 。

⑥ 如果滿足結(jié)束條件（誤差足夠好或達(dá)到最大循環(huán)次數(shù)）則退出，否則返回②。

在MATLAB中編程實(shí)現(xiàn)的基本粒子群算法基本函數(shù)為PSO，其調(diào)用格式如下：

[xm, fv] = PSO(fitness, N, c1, c2, w, M, D)

其中，fitness為待優(yōu)化的目標(biāo)函數(shù)，也稱適應(yīng)度函數(shù)。N是粒子數(shù)目，c1是學(xué)習(xí)因子1，c2是學(xué)習(xí)因子2，w是慣性權(quán)重，M是最大迭代次數(shù)，D是自變量的個數(shù)，xm是目標(biāo)函數(shù)取最小值時的自變量，fv是目標(biāo)函數(shù)的最小值。

使用MATLAB實(shí)現(xiàn)基本粒子群算法代碼如下：

function [xm, fv] = PSO(fitness, N, c1, c2, w, M, D)

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

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

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

% w慣性權(quán)重

% M最大迭代次數(shù)

% D搜索空間維數(shù)

% N初始化群體個數(shù)數(shù)目

%%%%%% 初始化種群的個體（可以在這里限定位置和速度的范圍） %%%%%%

format long;

for i = 1 : N

for j = 1 : D

    x(i, j) = randn;        % 隨機(jī)初始化位置

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

end

end

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

for i = 1 : N

p(i) = fitness(x(i, :));

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

end

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

for i = 1 : (N - 1)

if fitness(x(i, :)) < fitness(pg)

    pg = x(i, :);

end

end

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

for t = 1 : M

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 fitness(x(i, :)) < p(i)

        p(i) = fitness(x(i, :));

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

    end

    if p(i) < fitness(pg)

        pg = y(i, :);

    end

end

Pbest(t) = fitness(pg);

end

%%%%%% 最終給出計算結(jié)果 %%%%%%

disp('*************************************************')

disp('目標(biāo)函數(shù)取最小值時的自變量：')

xm = pg'

disp('目標(biāo)函數(shù)的最小值為：')

fv = fitness(pg)

disp('*************************************************')

將上面的函數(shù)保存到MATLAB可搜索路徑中，即可調(diào)用該函數(shù)。再定義不同的目標(biāo)函數(shù)fitness和其他輸入量，就可以用粒子群算法求解不同問題。

粒子群算法使用的函數(shù)有很多個，下面介紹兩個常用的適應(yīng)度函數(shù)。

1．Griewank函數(shù)

Griewank函數(shù)的MATLAB代碼如下：

function y = Griewank(x) % Griewank函數(shù)

% 輸入x，給定相應(yīng)的y值，在x = (0, 0, ……, 0)處有全局極小點(diǎn)0

[row, col] = size(x);

if row > 1

error('輸入的參數(shù)錯誤');

end

y1 = 1 / 4000 * sum(x .^ 2);

y2 = 1;

for h = 1 : col

y2 = y2 * cos(x(h) / sqrt(h));

end

y = y1 - y2 + 1;

y = - y;

繪制以上函數(shù)圖像的MATLAB代碼如下：

function DrawGriewank() % 繪制Griewank函數(shù)圖像

x = [-8 : 0.1 : 8];

y = x;

[X, Y] = meshgrid(x, y);

[row, col] = size(X);

for l = 1 : col

for h = l : row

    z(h, l) = Griewank([X(h, l), Y(h, l)]);

end

end

surf(X, Y, z);

shading interp

將以上代碼保存為DrawGriewank.m文件，并運(yùn)行上述代碼，得到Griewank函數(shù)圖像，如圖10-3所示。

圖10-3 Griewank函數(shù)圖像

2．Rastrigin函數(shù)

Rastrigin函數(shù)的MATLAB代碼如下：

function y = Rastrigin(x) % Rastrigin函數(shù)

% 輸入x，給定相應(yīng)的y值，在x = (0, 0, ……, 0)處有全局極小點(diǎn)0

[row, col] = size(x);

if row > 1

error('輸入的參數(shù)錯誤');

end

y = sum(x .^ 2 - 10 * cos(2 * pi * x) + 10);

y = - y;

繪制以上函數(shù)圖像的MATLAB代碼如下：

function DrawRastrigin()

x = [-4 : 0.05 : 4];

y = x;

[X, Y] = meshgrid(x, y);

[row, col] = size(X);

for l = 1 : col

for h = 1 : row

    z(h, l) = Rastrigin([X(h, l), Y(h, l)]);

end

end

surf(X, Y, z);

shading interp

將以上代碼保存為DrawRastrigin.m文件，并運(yùn)行上述代碼，得到Rastrigin函數(shù)圖像，如圖10-4所示。

圖10-4 Rastrigin函數(shù)圖像

例10-1：利用上文介紹的基本粒子群算法求解下列函數(shù)的最小值。

利用基本粒子群算法求解最小值，首先需要確認(rèn)不同迭代步數(shù)對結(jié)果的影響。設(shè)定題中函數(shù)的最小點(diǎn)均為0，粒子群規(guī)模為50，慣性權(quán)重為0.5，學(xué)習(xí)因子c1為1.5，學(xué)習(xí)因子c2為2.5，迭代步數(shù)分別取100、1000、10000。

在MATLAB中建立目標(biāo)函數(shù)代碼，并保存為fitness.m文件：

function F = fitness(x)

F = 0;

for i = 1 : 30

F = F + x(i)^2 + x(i) - 6

end

在MATLAB命令行窗口中依次輸入：

x = zeros(1, 30);

[xm1, fv1] = PSO(@fitness, 50, 1.5, 2.5, 0.5, 100, 30);

[xm2, fv2] = PSO(@fitness, 50, 1.5, 2.5, 0.5, 1000, 30);

[xm3, fv3] = PSO(@fitness, 50, 1.5, 2.5, 0.5, 10000, 30);

運(yùn)行以上代碼，比較目標(biāo)函數(shù)取最小值時的自變量，如表10-1所示。

表10-1 比較不同迭代步數(shù)下的目標(biāo)函數(shù)值和最小值

從表10-1中可以看出，迭代步數(shù)不一定與獲得解的精度成正比，即迭代步數(shù)越大，獲得解的精度不一定越高。這是因?yàn)榱Ｗ尤核惴ㄊ且环N隨機(jī)算法，同樣的參數(shù)也會算出不同的結(jié)果。

在上述參數(shù)的基礎(chǔ)上，保持慣性權(quán)重為0.5、學(xué)習(xí)因子c1為1.5、學(xué)習(xí)因子c2為2.5、迭代步數(shù)為100不變，粒子群規(guī)模分別取10、100和500，運(yùn)行以下MATLAB代碼：

x = zeros(1, 30);

[xm1, fv1] = PSO(@fitness, 10, 1.5, 2.5, 0.5, 100, 30);

[xm2, fv2] = PSO(@fitness, 100, 1.5, 2.5, 0.5, 100, 30);

[xm3, fv3] = PSO(@fitness, 500, 1.5, 2.5, 0.5, 100, 30);

比較目標(biāo)函數(shù)取最小值時的自變量，如表10-2所示。

表10-2 比較不同粒子群規(guī)模下的目標(biāo)函數(shù)值和最小值

從表10-2中可以看出，粒子群規(guī)模越大，獲得解的精度不一定越高。

綜合以上不同迭代步數(shù)和不同粒子群規(guī)模運(yùn)算得到的結(jié)果可知，在粒子群算法中，要想獲得精度高的解，關(guān)鍵是各個參數(shù)之間的匹配。