基于粒子群算法的分组背包MATLAB实现

抽空看了一段时间的粒子群算法，这里仅针对其应用于动态规划中的背包问题的情况做下总结归纳，其他应用可以之后想到了再添加。

如果这篇文章看完了的话，可以举一反三，将下面的代码稍微改动一下做一个二维函数的求最小值问题，点击链接可以直达哦。

一：分组背包问题简介

假设有3个组，每组有2个物品，每种物品有3种属性，价值、体积和重量。我们只有1个背包，从每组中选择1个物品（可以不选的情况第三章讨论）装入背包中，如何选择才能使背包中的物品总价值最大、总体积最小、且不超过规定重量呢？

由于上诉情况数较少，我们穷举就能列完：

假设背包总重量的限制为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 粒子运动计算

粒子速度的计算算是粒子群算法的精华了，其中有根据惯性权值计算粒子速度的公式：