基於粒子羣算法的分組揹包MATLAB實現

抽空看了一段時間的粒子羣算法，這裏僅針對其應用於動態規劃中的揹包問題的狀況作下總結概括，其餘應用能夠以後想到了再添加。

一：分組揹包問題簡介

假設有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的話，從窮舉中咱們能夠知道此時第一組物品二、第二組物品二、第三組物品1這樣選有總價值八、整體積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 粒子運動計算

粒子速度的計算算是粒子羣算法的精華了，其中有根據慣性權值計算粒子速度的公式：

\[v^{i+1} = wv^i + c1r1(p_{local}^i - x^i) + c2r2(p_{global} - x^i) \]

其中\(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.