數學建模方法-粒子群算法
一、引言
哈嘍大家好，有一段時間沒更新Blog了，最近身體不太舒服哈，今天開始繼續更了。言歸正傳，這次要講的是“粒子群算法”。這個算法是由兩個科學家在1995年，根據對鳥類捕食行為的研究所得到啟發而想出來的。好的，接下來讓我們開始吧。
二、鳥類捕食行為
鳥媽媽有7個鳥寶寶，有一天，鳥媽媽讓鳥寶寶們自己去找蟲子吃。於是鳥寶寶們開始了大範圍的捕食行為。一開始鳥寶寶們不知道哪裏可以找得到蟲子，於是每個鳥寶寶都朝著不同的方向獨自尋找。
但是為了能夠更快的找到蟲子吃，鳥寶寶們協商好，誰發現了蟲子，就互相說一聲。
找了一會，終於有一個鳥寶寶（稱之為小藍），似乎發現在他附近不遠處有蟲子的蹤跡。於是它傳話給其他鳥寶寶，其他鳥寶寶，收到消息後，邊開始改變軌跡，飛到小藍這邊。最終，隨著小藍越來越接近蟲子。其他蟲寶寶也差不多都聚集到了小藍這邊。最終，大家都吃到了蟲子。
三、粒子群算法
3.1 故事分析
鳥寶寶捕食的故事，正是這個粒子群算法存在的原因。因此，如果想更好的了解粒子群算法，我們就要來分析鳥寶寶捕食的故事。首先，我們來分析分析鳥寶寶們的運動狀態，即鳥寶寶自身是怎麽決定自己的飛翔速度和位置的。
(1) 吶首先，我們知道物體是具有慣性的，鳥寶寶在一開始飛翔的時候，無論它下一次想怎麽飛，往哪個方向飛，它都有一個慣性，它必須根據當前的速度和方向來進行下一步的調整。對吧，這個可以理解吧，因此，“慣性”——當前的速度$v_{current}$是一個因素。
(2) 其次，由於鳥寶寶長期捕食，因此鳥寶寶有經驗，它雖然不知道具體哪裏是有蟲子的存在，但是它能大概知道蟲子分布在哪裏。比如當鳥寶寶飛到貧瘠的地方，它肯定知道這裏是不會有蟲子的，因此，在鳥寶寶的心中，它每次飛，都會根據自己的經驗來找，比如以往蟲子分布的地區，它肯定優先對這部分的地方進行搜索。因此，自己的“認知
(3) 最後，每個鳥寶寶發現自己離蟲子更接近的時候，便會發出信號給同伴，在遇到這樣的信號，其余還在找的鳥寶寶們就會想著，同伴的位置更接近蟲子，我要往它那邊過去看看。我們稱之為“社會共享”，這也是鳥寶寶在移動時考慮的一個因素。
綜上所述，鳥寶寶每次在決定下一步移動的速度和方向時，腦子裏是由這三個因素影響的。我們想，如果能夠用一條公式來描述著三個因素的影響的話，那不就能寫出每個鳥寶寶的移動方向和速度麽。但是，每一個鳥寶寶，對這三個因素的考慮都是不一致的，比如對於捕食經驗不高的鳥寶寶，移動的時候會更看重同伴分享的信息，而對於捕食經驗高的鳥寶寶，則更看重自己的經驗。因此，我們的公式，在“認知”和“社會共享”這部分，是具有隨機性的。
3.2 公式表達
我們的粒子群算法是這樣的：
(1) 每一個鳥寶寶，都是我們的解，稱之為“粒子”，而我們的蟲子，就是我們問題的最優解。也就是說，鳥寶寶捕食過程，也就是所有的粒子在解空間內尋找到最優解的過程。
(2) 每一個粒子，都由一個fitness function來確定fitness value，以此來確定粒子位置的好壞。（這個其實就是模仿鳥寶寶的“經驗”判斷部分，通過fitness function來確定這個位置是不是所謂的貧瘠地或是蟲子可能出現的位置）。
(3) 每一個粒子被賦予了記憶功能，能夠記憶所搜尋過的最佳位置。
(4) 每一個粒子都有一個速度以及決定飛行的距離和方向，這個根據它本身的飛行經驗和同伴共享的信息所決定。
現在，在d維空間中，有n個粒子。其中：
粒子i的位置：
$X_{i} = (x_{i1}, x_{i2}, \cdot \cdot \cdot, x_{id})$
粒子i的速度：
$V_{i} = (v_{i1}, v_{i2}, \cdot \cdot \cdot, v_{id})$
粒子i所搜尋到的最好的位置（personal best）：
$Pbest_{id} = (pbest_{i1}, pbest_{i2}, \cdot \cdot \cdot, pbest_{id})$
種群所經歷的最好的位置(global best)：
$Gbest_{d} = (gbest_{1}, gbest_{2}, \cdot \cdot \cdot, gbest_{d})$
另外，我們要給我們的粒子速度和位置做一個限制，畢竟我們不希望鳥寶寶的速度過快或者以及飛出我們要搜尋的範圍。因此對於每個粒子i，有：
$X_{i} \in [X_{min}, X_{max}]$
$V_{i} \in [V_{min}, V_{max}]$
根據前面的分析，我們可以寫出下面兩條公式：
- 在d維空間中，粒子i的速度更新公式為：
$V_{i}^{k} = wV_{i}^{k-1}+c_{1}rand()\left (Pbest_{i} - X_{i}^{k-1} \right ) + c_{2}Rand() \left( Gbest_{i} - X_{i}^{k-1} \right )$
- 在d維空間中，粒子i的位置更新公式為：
$X_{i}^{k} = X_{id}^{k-1}+ V_{i}^{k-1}$
在上式中，上標k-1和k表示粒子從k-1次飛行操作到下一次飛行操作的過程。
(1) 我們先來看看速度更新公式，可以看出包含三部分，分別是我們前面分析的“慣性”、“認知”、“社會共享”這三大塊。而rand()和Rand()分別給“認知”和“社會共享”提供隨機性，即每個粒子對這兩部分的看重是不一樣的。其中c_{1}和c_{2}是我們自己加的，表示我們自己對這兩部分的分量的控制。另外w是一個慣性的權重，是用於調節對解空間的搜索範圍。關於這個w的出現還有一段歷史，大家感興趣的可以去查查論文，這裏就不細講了。
(2) 接著是位置更新公式，這部分很簡單，我們都知道$x = x_{0} + vt$。可是這裏為什麽少了個$t$呢，這裏我們可以簡單理解為從k-1次飛行到下一次飛行，耗費了一個單位時間$t$，因此就沒有$t$出現了。
好了，上面那兩個公式就是粒子群算法的核心了。
3.3 算法流程
(1) Initial:
初始化粒子群體（群體規模為n），每個粒子的信息包括隨機位置和隨機速度。
(2) Evaluation
根據fitness function，算出每個粒子的fitness value。
(3) Find the Pbest
對於每個粒子，將其當前的fitness value與歷史最佳的位置（Pbest）所對應的fitness value做比較。若當前的fitness value更高，則將當前的位置更新Pbest。
(4) Find the Gbest
對於每個粒子，將其當前的fitness value與群體歷史最佳的位置（Gbest）所對應的fitness value做比較。若當前的fitness value更高，則將當前的位置更新Gbest。
(5) Update the Velocity and Position:
根據公式更新每個粒子的速度和位置
(6) 如果未找到最佳值，則返回步驟2。（若達到了最佳叠代數量或者最佳fitness value的增量小於給定的閾值，則停止算法）
四、粒子群算法Matlab代碼示例
4.1 利用粒子群算法計算一元函數的最大值
%% 題目1：利用粒子群算法計算一元函數的最大值
%% I. 清空環境
clc
clear all
%% II. 繪制目標函數曲線圖
x = 1: 0.01: 2;
y = sin(10*pi*x) ./ x;
figure
plot(x, y)
hold on
%% III. 參數初始化
c1 = 1.49445;
c2 = 1.49445;
maxgen = 50; % 進化次數(叠代次數)
sizepop = 10; % 種群規模(粒子數數目)
dimension = 1; % 這裏因為是一元函數的求解，即一維，故列數為1
% 速度的邊界
Vmax = 0.5;
Vmin = -0.5;
% 種群的邊界
popmax = 2;
popmin = 1;
% 用於計算慣性權重，經驗值
ws = 0.9;
we = 0.4;
% 給矩陣預分配內存
pop = zeros(sizepop, dimension);
V = zeros(sizepop, dimension);
fitness = zeros(sizepop, 1);
yy = zeros(maxgen);
w = zeros(maxgen);
%% IV. 產生初始粒子和速度
for i = 1: sizepop
% 隨機產生一個種群
pop(i, :) = (rands(1) + 1) / 2 + 1; % 初始種群
V(i, :) = 0.5 * rands(1); % 初始化速度
% 計算適應度
fitness(i) = fun(pop(i, :));
end
%% V. 初始化Personal best和Global best
[bestfitness, bestindex] = max(fitness);
gbest = pop(bestindex,:); % Global best
pbest = pop; % 一開始的個體數據都是最佳的
fitnesspbest = fitness; % 個體最佳適應度值
fitnessgbest = bestfitness; % 全局最佳適應度值
%% VI. 叠代尋優
for i = 1: maxgen % 叠代循環
w(i) = ws - (ws - we) * (i / maxgen); % 讓慣性權重隨著叠代次數而動態改變，控制搜索精度
for j = 1: sizepop % 遍歷所有粒子
% 速度更新
V(j, :) = w(i)*V(j, :) + c1*rand*(pbest(j, :) - pop(j, :)) + c2*rand*(gbest - pop(j, :));
if V(j, :) > Vmax
V(j, :) = Vmax;
end
if V(j, :) < Vmin
V(j, :) = Vmin;
end
% 種群更新(位置更新)
pop(j, :) = pop(j, :) + V(j, :);
if pop(j, :) > popmax
pop(j, :) = popmax;
end
if pop(j, :) < popmin
pop(j, :) = popmin;
end
% 適應度值更新
fitness(j) = fun(pop(j, :));
end
for j = 1: sizepop
% 個體最優更新
if fitness(j) > fitnesspbest(j)
pbest(j, :) = pop(j, :);
fitnesspbest(j) = fitness(j);
end
%群體最優更新
if fitness(j) > fitnessgbest
gbest = pop(j, :);
fitnessgbest = fitness(j);
end
end
yy(i) = fitnessgbest; % 記錄每次叠代完畢的群體最優解
end
%% VII. 輸出結果並繪圖
[fitnessgbest gbest] % 輸出數據
figure
plot(yy)
title(‘最優個體適應度‘, ‘fontsize‘,12);
xlabel(‘進化代數‘, ‘fontsize‘, 12);
ylabel(‘適應度‘, ‘fontsize‘,12);
其中，適應值函數被封裝到fun.m中，如下：
function y = fun(x) % 函數用於計算粒子適應度值 %x input 輸入粒子 %y output 粒子適應度值 y = sin(10 * pi * x) / x;
運行上述代碼，得到的數據和圖如下：
ans =
0.9528 1.0490
可以看到，紅點處正是我們函數的最大值處。
4.2 利用粒子群算法計算二元函數的最大值
%% 題目2： 利用粒子群算法計算二元函數的最大值
%% I. 清空環境
clc
clear
%% II. 繪制目標函數曲線
figure
[x, y] = meshgrid(-5: 0.1: 5, -5: 0.1: 5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20;
mesh(x,y,z)
hold on
%% III. 參數初始化
c1 = 1.49445;
c2 = 1.49445;
maxgen = 1000; % 進化次數
sizepop = 100; % 種群規模
dimension = 2; % 這裏因為是二元函數的求解，即二維，故列數為2
% 速度的邊界
Vmax = 1;
Vmin = -1;
% 種群的邊界
popmax = 5;
popmin = -5;
% 用於計算慣性權重，經驗值
ws = 0.9;
we = 0.4;
% 給矩陣預分配內存
pop = zeros(sizepop, dimension);
V = zeros(sizepop, dimension);
fitness = zeros(sizepop, 1);
yy = zeros(maxgen);
w = zeros(maxgen);
%% IV. 產生初始粒子和速度
for i = 1: sizepop
% 隨機產生一個種群
pop(i, :) = 5 * rands(1, 2); % 初始種群
V(i, :) = rands(1, 2); % 初始化速度
% 計算適應度
fitness(i) = fun(pop(i, :));
end
%% V. 初始化Personal best和Global best
[bestfitness, bestindex] = max(fitness);
gbest = pop(bestindex, :); % Global best
pbest = pop; % 個體最佳
fitnesspbest = fitness; % 個體最佳適應度值
fitnessgbest = bestfitness; % 全局最佳適應度值
%% VI. 叠代尋優
for i = 1: maxgen
w(i) = ws - (ws - we) * (i / maxgen); % 讓慣性權重隨著叠代次數而動態改變，控制搜索精度
for j = 1: sizepop
% 速度更新
V(j, :) = w(i)*V(j, :) + c1*rand*(pbest(j, :) - pop(j, :)) + c2*rand*(gbest - pop(j, :));
for k = 1: dimension
if V(j, k) > Vmax
V(j, k) = Vmax;
end
if V(j, k) < Vmin
V(j, k) = Vmin;
end
end
% 種群更新(位置更新)
pop(j, :) = pop(j, :) + V(j, :);
for k = 1: dimension
if pop(j, k) > popmax
pop(j, k) = popmax;
end
if pop(j, k) < popmin
pop(j, k) = popmin;
end
end
% 適應度值更新
fitness(j) = fun(pop(j, :));
end
for j = 1: sizepop
% 個體最優更新
if fitness(j) > fitnesspbest(j)
pbest(j, :) = pop(j, :);
fitnesspbest(j) = fitness(j);
end
%群體最優更新
if fitness(j) > fitnessgbest
gbest = pop(j, :);
fitnessgbest = fitness(j);
end
end
yy(i) = fitnessgbest; % 記錄每次叠代完畢的群體最優解
end
%% VII.輸出結果
[fitnessgbest, gbest]
plot3(gbest(1), gbest(2), fitnessgbest, ‘bo‘,‘linewidth‘, 1.5)
figure
plot(yy)
title(‘最優個體適應度‘, ‘fontsize‘, 12);
xlabel(‘進化代數‘, ‘fontsize‘, 12);
ylabel(‘適應度‘, ‘fontsize‘, 12);
其中，適應值函數被封裝到fun.m中，如下：
function y = fun(x) %函數用於計算粒子適應度值 %x input 輸入粒子 %y output 粒子適應度值 y = x(1).^2 + x(2).^2 - 10*cos(2*pi*x(1)) - 10*cos(2*pi*x(2)) + 20;
運行上述代碼，得到的數據和圖如下：
ans = 80.7066 4.5230 4.5230
可以看到，圖中標註的地方是我們求得的最大值處。其實我們知道，對於這個函數，因為是對稱的，所以4個點都是同樣的最大值，這就是粒子群算法的缺點了。很容易陷入局部最優解。因為我們的粒子群算法其實並不知道什麽是最優解，它只是讓粒子不斷的找一個相對之前所有的解都是最好的一個解。所以，這樣的粒子群算法是有局限性的，用的時候要根據場合來選擇。
數學建模方法-粒子群算法