抽空看了一段時間的粒子群演算法，這裡僅針對其應用於動態規劃中的揹包問題的情況做下總結歸納，其他應用可以之後想到了再新增。
一：分組揹包問題簡介
假設有3個組，每組有2個物品，每種物品有3種屬性，價值、體積和重量。我們只有1個揹包，從每組中選擇1個物品（可以不選的情況第三章討論）裝入揹包中，如何選擇才能使揹包中的物品總價值最大、總體積最小、且不超過規定重量呢？
| 物品/分組 | 第一組 | 第二組 | 第三組 |
|---|---|---|---|
| 物品1價值 | 1 | 2 | 3 |
| 物品2價值 | 3 | 2 | 1 |
| 物品/分組 | 第一組 | 第二組 | 第三組 |
|---|---|---|---|
| 物品1體積 | 1 | 2 | 1 |
| 物品2體積 | 2 | 1 | 2 |
| 物品/分組 | 第一組 | 第二組 | 第三組 |
|---|---|---|---|
| 物品1質量 | 20 | 50 | 40 |
| 物品2質量 | 30 | 40 | 20 |
由於上訴情況數較少，我們窮舉就能列完：
| 每組選第幾物品 | 揹包總價值 | 揹包總體積 | 揹包總重量 |
|---|---|---|---|
| 1 1 1 | 6 | 4 | 110 |
| 1 1 2 | 4 | 5 | 90 |
| 1 2 1 | 6 | 3 | 100 |
| 1 2 2 | 4 | 4 | 80 |
| 2 1 1 | 8 | 5 | 120 |
| 2 1 2 | 6 | 6 | 100 |
| 2 2 1 | 8 | 4 | 110 |
| 2 2 2 | 6 | 5 | 90 |
假設揹包總重量的限制為110的話，從窮舉中我們可以知道此時第一組物品2、第二組物品2、第三組物品1這樣選有總價值8、總體積4這樣的最優解。下面將介紹如何用粒子群演算法求解此類問題。
二：粒子群演算法解分組揹包問題
一般說到揹包問題都會想到用動態規劃DP的方法去解決，DP確實很精煉簡潔，但是狀態轉移方程的理解不易。粒子群演算法過程長了些，但其思想很樸素，更自然。
先說說學習理解了演算法過後，對這演算法運算過程的感性描述：隨機產生了一堆粒子，每個粒子代表揹包的一種情況（選了哪3個物品），初始粒子全都是區域性最優粒子，然後從這堆粒子中按優劣選出非劣粒子。之後進入迭代過程，每次迭代中，從非劣粒子中隨機選一個作為全域性最優粒子，按照公式計算粒子速度（跟較優粒子和全域性最優粒子息息相關），使每個粒子產生移動（也就是有概率替換代表的物品），如果移動後的粒子比之前的粒子更優則替代原來的區域性最優粒子，然後把非劣粒子和區域性最優粒子放一起，再按優劣選出非劣粒子，去除重複的非劣粒子，進入下一次迭代。
注意加粗的部分都是演算法的關鍵步驟！
自己設定迭代次數，粒子群演算法收斂很快，最後如果有最優解的話，我們會得到一堆非劣解粒子，選擇裡面價值最大且體積小的情況就可以啦！
重點是看程式是如何實現的，下面將完整的MATLAB程式碼分解成2.1~2.7七部分講解。
2.1 輸入引數、固定引數初始化
clear, clc, close all;
%% 輸入引數、固定引數初始化
P = [1 2 3; 3 2 1]; % 各組物品價值
V = [1 2 1; 2 1 2]; % 各組物品體積
M = [20 50 40; 30 40 20]; % 各組物品重量
group = size(P, 2); % 組數
nitem = size(P, 1); % 每組物品數
weight = 120; % 揹包最大重量限制
xsize = 50; % 粒子數
ITER = 200; % 迭代次數
c1 = 0.8; c2 = 0.8; % 常數
wmax = 1.2; wmin = 0.1; % 慣性權重相關常數
v = zeros(xsize, group); % 粒子速度初始化
2.2 粒子群位置、適應度、最佳位置、最佳適應度初始化
隨機產生粒子群\(x\)，表示每組選哪個物品（其實就是產生了一群揹包），比如\(x_i = [1, 2, 1]\)的時候，表示第\(i\)個粒子取第一組第1個物品、第二組第1個物品、第三組第1個物品
注意粒子與粒子群位置是一個意思。有時候可能會混用。
適應度其實就是用粒子群位置算出其代表的揹包的價值、體積、質量。
%% 粒子群位置、適應度、最佳位置、最佳適應度初始化
x = randi(nitem, xsize, group); % 隨機粒子群位置（表示每組選哪個物品）
% 粒子群適應度
xp = zeros(1, xsize); % 粒子群價值
xv = zeros(1, xsize); % 粒子群體積
xm = zeros(1, xsize); % 粒子群重量
for i = 1 : xsize
for j = 1 : group
xp(i) = xp(i) + P(x(i, j), j);
xv(i) = xv(i) + V(x(i, j), j);
xm(i) = xm(i) + M(x(i, j), j);
end
end
bestx = x; % 粒子群位置最佳值
bestp = xp; bestv = xv; bestm = xm; % 粒子群最佳適應度
2.3 初始篩選非劣解
第一次篩選非劣解，之後每次迭代都會重新篩選一次。其中的判斷條件很重要，可以根據問題的限制而改變。這裡就是判斷每個粒子是否比別的所有粒子都更符合要求（價值大且體積小）。
如果還限制了揹包體積大小的話，還需要在判斷中加上對粒子體積的限制。這裡只限制了揹包重量。
%% 初始篩選非劣解
cnt = 1;
for i = 1 : xsize
fljflag = 1;
for j = 1 : xsize
if j ~= i
if (xp(j) > xp(i) && xv(j) < xv(i)) ||... % i粒子劣解
(xp(j) > xp(i) && xv(j) == xv(i)) ||...
(xp(j) == xp(i) && xv(j) < xv(i)) || (xm(i) > weight)
fljflag = 0;
break;
end
end
end
if fljflag == 1
flj(cnt, :) = [xp(i), xv(i), xm(i)]; % 非劣解適應度
fljx(cnt, :) = x(i, :); % 非劣解位置
cnt = cnt + 1;
end
end
for niter = 1 : ITER % 迭代開始，粒子開始運動
rnd = randi(size(flj, 1), 1, 1);
gbest = fljx(rnd, :); % 粒子全域性最優解
2.4 粒子運動計算
粒子速度的計算算是粒子群演算法的精華了，其中有根據慣性權值計算粒子速度的公式：
其中\(w\)為慣性權值，\(c1\)和\(c2\)為常數，\(r1\)和\(r2\)為[0，1]間服從均勻分佈的隨機數，\(p_{local}^i\)是區域性最優粒子，\(p_{global}\)是全域性最優粒子，注意只有一個所以沒有上標\(i\)。
慣性權值的取值跟迭代次數有關，這裡我們採用\(w = wmax-(wmax - wmin) * niter / iterall\)這樣的計算方法。相關慣性權值的計算也是粒子群演算法研究的熱點，慣性權值變化大，粒子速度快，位置變換也快，慣性權值取得好的話可以使粒子群更快的收斂到全域性最優！
在用速度對每個粒子位置進行更新時，注意兩個問題，一是粒子位置可能取到大於物品總數的值了，所以這裡我們取了一個模值；二是粒子運動後可能會產生負的值，總不能取到第-2個物品吧，所以在這裡的處理方法是對負數取一個小於等於物品數的隨機正整數。最後記得對運動後的粒子位置值向上取整。
%% 粒子運動計算
w = wmax - (wmax - wmin) * niter / ITER; % 慣性權值更新
r1 = rand(1,1); r2 = rand(1,1); % 產生[0, 1]間均勻分佈隨機值
for i = 1 : xsize
v(i, :) = w * v(i, :) + c1 * r1 * (bestx(i, :) - x(i, :)) + c2 * r2 * (gbest - x(i, :)); % 粒子速度
x(i, :) = x(i, :) + v(i, :);
x(i, :) = rem(x(i, :), nitem); % 不能超過物品數
index = find(x(i, :) <= 0); % 運動後小於零了，用新的隨機整數代替
if ~isempty(index)
x(i, index) = randi(nitem, 1, length(index));
end
x(i, :) = ceil(x(i, :));
end
2.5 當前粒子群適應度、最佳位置、最佳適應度
粒子經過了運動到了新的位置，當然要把新粒子和舊粒子拿出來比一比，如果新粒子的適應度比舊粒子好的話，就更新區域性最佳粒子位置咯。這裡需要考慮一個問題，新粒子好的話更新，不好的話不更新，但是遇到新舊各有所長時怎麼辦呢？其實就相當於有兩個揹包，一個價值更高，另一個體積更小，那到底取哪一個呢？那就隨機取一個吧。
%% 當前粒子群適應度、最佳位置、最佳適應度
xp_cur = zeros(1, xsize);
xv_cur = zeros(1, xsize);
xm_cur = zeros(1, xsize);
for i = 1 : xsize
for j = 1 : group
xp_cur(i) = xp_cur(i) + P(x(i, j), j);
xv_cur(i) = xv_cur(i) + V(x(i, j), j);
xm_cur(i) = xm_cur(i) + M(x(i, j), j);
end
end
for i = 1 : xsize
if (xp_cur(i) > xp(i) && xv_cur(i) < xv(i)) ||... % 如果當前粒子適應度更好
(xp_cur(i) > xp(i) && xv_cur(i) == xv(i)) ||...
(xp_cur(i) == xp(i) && xv_cur(i) < xv(i)) || (xm(i) > weight)
bestx(i, :) = x(i, :); % 粒子最佳位置更新
bestp(i) = xp_cur(i); bestv(i) = xv_cur(i); bestm(i) = xm_cur(i); % 粒子最佳適應度更新
elseif (xp_cur(i) < xp(i) && xv_cur(i) > xv(i)) ||... % 如果原始粒子適應度更好，不做處理
(xp_cur(i) < xp(i) && xv_cur(i) == xv(i)) ||...
(xp_cur(i) == xp(i) && xv_cur(i) > xv(i)) || (xm_cur(i) > weight)
continue;
else % 當前與原始粒子適應度各有所長，隨機更新
if rand(1,1) < 0.5
bestx(i, :) = x(i, :); % 粒子最佳位置更新
bestp(i) = xp_cur(i); bestv(i) = xv_cur(i); bestm(i) = xm_cur(i); % 粒子最佳適應度更新
end
end
end
xp = xp_cur; xv = xv_cur; xm = xm_cur; % 新粒子變成舊粒子
2.6 粒子群最佳位置、最佳適應度合併後再篩選非劣解
與第一次做非劣解篩選的步驟基本一樣，只是加了一個合併操作，把區域性最佳粒子與非劣解粒子合併在一起，然後再篩選一波非劣解粒子。
%% 粒子群最佳位置、最佳適應度合併後再篩選非劣解
bestxmerge = [bestx; fljx];
pmerge = [bestp, flj(:, 1)'];
vmerge = [bestv, flj(:, 2)'];
mmerge = [bestm, flj(:, 3)'];
n = size(flj, 1);
flj = [];
fljx = [];
cnt = 1;
for i = 1 : xsize + n
fljflag = 1;
for j = 1 : xsize + n
if j ~= i
if (pmerge(j) > pmerge(i) && vmerge(j) < vmerge(i)) ||...
(pmerge(j) > pmerge(i) && vmerge(j) == vmerge(i)) ||...
(pmerge(j) == pmerge(i) && vmerge(j) < vmerge(i)) || (mmerge(i) > weight)
fljflag = 0;
break;
end
end
end
if fljflag == 1
flj(cnt, :) = [pmerge(i), vmerge(i), mmerge(i)]; % 非劣解適應度
fljx(cnt, :) = bestxmerge(i, :); % 非劣解位置
cnt = cnt + 1;
end
end
2.7 去掉重複非劣解
這一步也很重要，實現方法也有很多，隨便選一種把重複的非劣解去掉就可以咯。
%% 去掉重複非劣解
issame = zeros(cnt - 1, 1);
for i = 1 : cnt - 1
for j = i + 1 : cnt - 1
if ~issame(j)
issame(j) = (length(find(fljx(j, :) == fljx(i, :))) == group);
end
end
end
flj(find(issame == 1), :) = [];
fljx(find(issame == 1), :) = [];
end
figure;
plot(flj(:, 1), flj(:, 2), 'ro');
title('粒子群結果'); xlabel('總價值'); ylabel('總體積');
2.8 結論分析
可以看到總價值8和總體積4的粒子（揹包）是最優揹包，符合我們最開始窮舉的結論。
三：粒子群演算法解其他揹包問題
3.1 分組揹包每組最多選一個物品
當求解分組揹包問題的另一種情況：每組物品最多選一個，也就是每組物品可選可不選時，我們可以把輸入引數矩陣做個小處理，新增1個物品，其價值、體積、質量都置為0就可以了。
| 物品/分組 | 第一組 | 第二組 | 第三組 |
|---|---|---|---|
| 物品1價值 | 1 | 2 | 3 |
| 物品2價值 | 3 | 2 | 1 |
| 物品3價值 | 0 | 0 | 0 |
| 物品/分組 | 第一組 | 第二組 | 第三組 |
|---|---|---|---|
| 物品1體積 | 1 | 2 | 1 |
| 物品2體積 | 2 | 1 | 2 |
| 物品3體積 | 0 | 0 | 0 |
| 物品/分組 | 第一組 | 第二組 | 第三組 |
|---|---|---|---|
| 物品1質量 | 20 | 50 | 40 |
| 物品2質量 | 30 | 40 | 20 |
| 物品3質量 | 0 | 0 | 0 |
四：常規DP演算法解分組揹包問題
4.1 C++ 常規DP解法
#include <iostream>
using namespace std;
int v[110][110],w[110][110],s[110];
int f[110];
int main()
{
int N,V;
cin >> N >> V;
for (int i = 1; i <= N; i ++)
{
cin >> s[i];
for (int j = 1; j <= s[i]; j ++)
cin >> v[i][j] >> w[i][j];
}
for (int i = 1; i <= N; i ++)
for (int j = V; j > 0; j --)
for (int k = 1; k <= s[i]; k ++)
if (j >= v[i][k])
f[j] = max(f[j],f[j - v[i][k]] + w[i][k]);
cout << f[V];
return 0;
}
/*
輸入：
3 5
2
1 2
2 4
1
3 4
1
4 5
輸出：
8
*/
4.2 MATLAB 粒子群解法
這道題目中少了一個重量的維度。我們只用把上面粒子群程式的輸入矩陣改為如下形式，然後把對重量的限制改為對體積的限制就好了。
P = [2 4 5; 4 0 0; 0 0 0]; % 各組物品價值
V = [1 3 4; 2 0 0; 0 0 0]; % 各組物品體積
weight = 5; % 揹包最大體積限制
可以看到總價值最高的粒子（揹包）也是8。
五：參考文獻
[1] 鬱磊,史峰,王輝,等.MATLAB智慧演算法30個案例分析(第2版)[M].北京:北京航空航天大學出版社,2015.