基于粒子群算法的多目标搜索算法（附代码案例）

一、理论基础

        在实际工程优化问题中,多数问题是多目标优化问题。相对于单目标优化问题,多目标优化问题的显著特点是优化各个目标使其同时达到综合的最优值。然而,由于多目标优化问题的各个目标之间的相互影响,多个目标难以同时达到最优,因此，一般适用于单目标问题的方法难以用于多目标问题的求解。
        多目标优化问题很早就引起了人们的重视，现已经发展出多种求解多目标优化问题的方法。多目标优化问题求解中最重要的概念是非劣解和非劣解集，两者的定义如下：
        非劣解(noninferior solution)：在多目标优化问题的可行域中存在一个问题解,若不存在另一个可行解,使得一个解中的目标全部劣于该解,则该解称为多目标优化问题的非劣解。所有非劣解的集合叫做非劣解集(noninferior set)。

        在求解实际问题中，过多的非劣解是无法直接应用的，决策者只能选择其中最满意的一个非劣解作为最终解。最终解选择主要有三种方法。第一种是求非劣解的生成法，包括加权法、约束法、加权法和约束法结合的混合法以及多目标遗传算法，即先求出大量的非劣解，构成非劣解的一个子集，然后按照决策者的意图找出最终解。第二种为交互法，主要为求解线性约束多目标优化的Geoffrion 法，不先求出很多的非劣解，而是通过分析者与决策者对话的方式，逐步求出最终解。第三种是事先要求决策者提供目标之间的相对重要程度，算法以此为依据，将多目标问题转化为单目标问题进行求解。
        目前，采用多目标进化算法求解多目标问题已成为进化计算领域中的一个热门方向，粒子群优化、蚁群算法、人工免疫系统、分布估计算法、协同进化算法、进化算法等一些新的进化算法陆续被用于求解多目标优化问题。本案例采用多目标粒子群算法求解多目标背包问题。

二、案例背景

1、问题描述

        假设存在五类物品，每类物品中又包含四种具体物品，现要求从这五种类别物品中分别选择一种物品放入背包中，使得背包内物品的总价值最大，总体积最小，并且背包的总质量不超过92 kg。背包问题的数学模型为：

        其中，Px 表示背包内物品价值；Rx 表示背包内物品体积；C 表示物品质量；X 为选择物品。P 为每个物品的价值，R 为每个物品的体积。P，R，C 的取值如下：

 2、算法流程

        基于粒子群算法的多目标搜索算法流程如下图所示。其中，种群初始化模块随机初始化粒子的位置 α 和速度 v，适应度值计算模块根据适应度值计算公式计算个体适应度值，粒子最优更新模块根据新的粒子位置更新个体最优粒子。非劣解集更新模块根据新粒子支配关系筛选非劣解。粒子速度和位置更新模块根据个体最优粒子位置和全局粒子位置更新粒子速度和位置。

         每个个体的适应度值有两个，即价值和体积，同时个体需满足质量约束。

3、筛选非劣解集

        筛选非劣解集主要分为初始筛选非劣解集和更新非劣解集。初始筛选非劣解集是指在粒子初始化后,当一个粒子不受其他粒子支配(即不存在其他粒子的P.,R均优于该粒子)时，把粒子放入非劣解集中,并且在粒子更新前从非劣解集中随机选择一个粒子作为群体最优粒子。更新非劣解集是指当新粒子不受其他粒子以及当前非劣解集中粒子的支配时,把新粒子放人非劣解集中,并且每次粒子更新前都从非劣解集中随机选择一个粒子作为群体最优粒子。

4、粒子速度和位置更新

        粒子更新公式如下：

        其中，w为惯性权重；r1 和 r2 为分布于 [ 0，1 ] 区间的随机数；k 是当前迭代次数；Pk(id) 为个体最优粒子位置；Pk(gd) 为全局最优粒子位置；c1 和 c2 为常数；V 为粒子速度；X 为粒子位置。

5、粒子最优

        粒子最优包括个体最优粒子和群体最优粒子，其中个体最优粒子的更新方式是从当前新粒子和个体最优粒子中选择支配粒子，当两个粒子都不是支配粒子时，从中随机选择一个粒子作为个体最优粒子。群体最优粒子为从非劣解集中随机选择的一个粒子。

三、MATLAB 程序实现

        根据多目标搜索算法原理，在 MATLAB 中实现基于粒子群算法的多目标搜索算法。

1、种群初始化

        初始化种群并且计算初始粒子的适应度值，程序代码如下：

Dim = 5;                            % 粒子维数
xSize = 50;                         % 种群个数

% 初始化粒子位置和速度
x = unidrnd(objnum, xSize, Dim);    % 粒子初始化 从具有标量最大值n的离散均匀分布生成随机数数组
% x = rand(xSize, Dim);
% x = ceil(x*objnum);               % 与unidrnd等价
v = zeros(xSize, Dim);              % 速度初始化

xbest = x;                          % 个体最佳值
gbest = x(1, :);                    % 粒子群最佳位置

% 粒子适应度值 
px = zeros(1, xSize);               % 粒子价值目标
rx = zeros(1, xSize);               % 粒子体积目标
cx = zeros(1, xSize);               % 重量约束

% 最优值初始化
pxbest = zeros(1, xSize);           % 粒子最优价值目标
rxbest = zeros(1, xSize);           % 粒子最优体积目标
cxbest = zeros(1, xSize);           % 记录重量，以求约束

% 计算初始种群适应度值
for i = 1: xSize
    for j = 1: Dim
        px(i) = px(i) + P(x(i, j), j);              % 粒子价值
        rx(i) = rx(i) + R(x(i, j), j);              % 粒子体积
        cx(i) = cx(i) + C(x(i, j), j);              % 粒子质量
    end
end

% 粒子最优位置
pxbest = px; rebest = rx; cxbest = cx;              % 粒子历史最优

2、种群更新

        种群更新根据全局最优粒子和个体最优粒子更新当前个体的速度和位置，其中全局最优粒子为非劣解集中随机选取的粒子。程序代码如下：

    % 从非劣解中选择粒子作为全局最优解
    s = size(fljx, 1);
    index = randi(s, 1, 1);         % randi返回一个由1到s的伪随机整数组成的1*1数组
    gbest = fljx(index, :);

    %% 群体更新
    for i=1:xSize
        % 速度更新
        v(i, :) = w*v(i,:) + c1*rand(1,1)*(xbest(i,:)-x(i,:)) + c2*rand(1,1)*(gbest-x(i,:);
        % 位置更新
        x(i, :) = x(i, :) + v(i, :);
        x(i, :) = rem(x(i, :), objnum) / double(objnum);    % rem 返回第一个数除以第二个数的余数
        indexl = fin(x(i, :) <= 0);
        if ~isempty(indexl)                                 % isempty 判断数组是否为空
            x(i, indexl) = rand(size(indexl));
        end
        x(i, :) = ceil(objnum * x(i, :));                   % ceil 朝正无穷大四舍五入
    end

3、更新个体最优粒子

        根据新粒子和当前最优粒子的支配关系，更新个体最优粒子，即当两个粒子存在支配粒子时，选择支配粒子，否则从中随机选取一个粒子作为新的个体最优粒子。程序代码如下：

    %% 更新粒子历史最佳
    for i=1:xSize
        %现在的支配原有的，替代原有的
        if ((px(i)<pxPrior(i)) && (rx(i)>rxPrior(i))) || ...
                    ((abs(px(i)-pxPrior(i))<tol) && (rx(i)>rxPrior(i))) || ...
                    ((px(i)<pxPrior(i)) && (abs(rx(i)-rxPrior(i))<tol)) || (cx(i)>weight)
                xbest(i, :) = x(i, :);
                pxbest(i)=pxPrior(i);rxbest(i)=rxPrior(i);cxbest(i)=cxPrior(i);
        end
        %彼此不受支配，随机决定
        if ~(((px(i)<pxPrior(i)) && (rx(i)>rxPrior(i))) || ...
                ((abs(px(i)-pxPrior(i))<tol) && (rx(i)>rxPrior(i))) || ...
                ((px(i)<pxPrior(i)) && (abs(rx(i)-rxPrior(i))<tol)) || (cx(i)>weight) ) ...
                &&  ~(((pxPrior(i)<px(i)) && (rxPrior(i)>rx(i))) || ...
                ((abs(pxPrior(i)-px(i))<tol) && (rxPrior(i)>rx(i))) || ...
                ((pxPrior(i)<px(i)) && (abs(rxPrior(i)-rx(i))<tol)) || (cxPrior(i)>weight))
            if rand(1, 1) < 0.5
                xbest(i, :) = x(i, :);
                pxbest(i)=pxPrior(i); rxbest(i)=rxPrior(i); cxbest(i)=cxPrior(i);
            end
        end
    end

4、非劣解筛选

        非劣解集筛选分为两步，第一步是把新非劣解集和旧非劣解集合并，得到新的非劣解集。程序代码如下：

    % 新个体非劣解
    px = pxPrior;
    rx = rxPrior;
    cx = cxPrior;
    % 更新升级非劣解集合
    s = size(flj, 1);                   % 目前非劣解集合中元素个数
    % 先将非劣解集合和xbest合并
    pppx = zeros(1, s + xSize);
    rrrx = zeros(1, s + xSize);
    cccx = zeros(1, s + xSize);
    
    pppx(1: xSize) = pxbest; pppx(xSize+1: end) = flj(:, 1)';
    rrrx(1: xSize) = pxbest; rrrx(xSize+1: end) = flj(:, 2)';
    cccx(1: xSize) = pxbest; cccx(xSize+1: end) = flj(:, 3)';
    
    xxbest = zeros(s + xSize, Dim);
    xxbest(1: xSize, :) = xbest;
    xxbest(xSize + 1: end, :) = fljx;

        第二步是根据非劣解集中的支配关系，筛选出新的非劣解集。程序代码如下：

    %筛选非劣解
    flj = [];
    fljx = [];
    k = 0;
    tol = 1e-7;
    for i = 1: xSize + s
        flag = 0;
        % 判断该点是否非劣
        for j = 1: xSize + s
            if j ~= i
                if ((pppx(i)<pppx(j)) && (rrrx(i)>rrrx(j))) || ...
                        ((abs(pppx(i)-pppx(j))<tol) && (rrrx(i)>rrrx(j))) || ...
                        ((pppx(i)<pppx(j)) && (abs(rrrx(i)-rrrx(j))<tol)) || ...
                        (cccx(i)>weight)        % 有一次被支配
                    flag = 1;
                    break;
                end
            end
        end
        % 判断有无被支配
        if flag == 0
            k = k + 1;
            % 记录非劣解
            flj(k, 1)=pppx(i); flj(k, 2)=rrrx(i); flj(k, 3)=cccx(i);
            % 记录非劣解位置
            fljx(k, :) = xxbest(i, :);
        end
    end

四、仿真结果

        从每类物品中选择一个物品放入背包，使背包的总价值最大，体积最小，并且背包总质量小于92 kg。粒子群算法参数为：粒子个数为50，迭代次数为200，最终得到的非劣解在目标空间中的分布如下图所示：

        由图知，算法搜索到的非劣解构成了 Pareto 面，算法搜索取得了很好的效果。

        多目标搜索算法相对于单目标算法来说，更加贴近于实际问题，求解结果更具有参考价值。通过多目标搜索算法最终得到的不是一个最优解，而是一个非劣解集，需要从非劣解集中根据实际问题的需要选择一个解作为该问题的最终解。常用的基于进化算法的多目标搜索算法除了本案例介绍的方法之外，还有基于遗传算法的多目标搜索算法，基于免疫算法的多目标搜索算法等。

代码文件和数据：

链接：https://pan.baidu.com/s/1YhwImXhtIUwIrZlYHT5jCg 
提取码：1ret