多目标粒子群算法(MOPSO)的原理和matlab实现

        算法原理部分参考文献基于改进多目标粒子群算法的配电网储能选址定容

0.前言

        初学者面对多目标优化问题可能比较困难，写下这篇博客记录一下自己学习的心得，希望能和大家一起交流学习。

        采用粒子群求单目标优化问题的原理很好理解，就是通过对粒子群的速度和位置不断来更新粒子群的最优适应度(也就是目标函数)，达到寻优的目的。但是多目标优化问题就比较难办了，由于目标函数有多个，如果同样使用粒子群算法，那么适应度怎么设置？怎么确定全体最优的粒子？这些问题都是比较棘手的。

        目前采用粒子群算法求解多目标问题主要有两大类处理方式：

        1.通过加权求和、灰色关联度分析、topsis、层次分析法等方法对多个目标函数进行处理，将多目标问题转换为单目标问题，然后仿照单目标问题进行求解；

        2.在迭代时以一定规则选择全体最优及个体最优的粒子，然后通过迭代求取pareto解集，最后根据需求从pareto解集中选取一个最合适的方案。

        第一种方式本质上还是单目标优化问题，这篇主要介绍一下第二种方法：

一、一些名词的解释

1.支配与非支配

        我们知道在粒子群算法中一个粒子就代表优化问题的一个解(也可以说是方案)，如果说在一个多目标优化问题中，解1得到的所有目标函数全部优于解2，那么就说解2被解1支配。啥意思呢，举个例子，比如咱们的优化目标是选房子，目标函数包括房子租金和面积。A租金1800/月，面积100m²，B租金2800/月，面积40m²，显然A完全优于B，那么就说解B被解A支配；再有C租金1600/月，面积80m²，C的租金低于A，但面积比A小，这时候就说不上来C和A哪个更优，也就是C是A(A是C)的非支配解。

2.pareto最优解与pareto最优解集

       清楚了支配与非支配的含义，那么pareto最优解(也叫非劣解)就不难理解了：如果说对于优化问题的一个解A，不存在其他的解可以支配解A，那么解A就是一个pareto最优解。显然pareto最优解不只一个，由这些解组成的集合就称为pareto最优解集。pareto最优解集一个重要的特点就是集合中所有解的互相之间都是非支配的关系，也就是解集中不存在任何一个解完全优于其他解(要和平共处)。所以，和PSO最后要求出唯一一个最优解不同，MOPSO的目标是求出一个互相之间为非支配关系的解的集合，也就是pareto最优解集。

3.外部存档

        因为粒子群算法在迭代过程中各个粒子的速度和位置都是不断变化的，适应度(目标函数)也随之变化，所以需要一般都需要采用一个外部存档将pareto最优解的数据存储下来。

二、算法原理和步骤

step1 种群初始化

        根据所求问题的具体情况，初始化所有粒子的速度和位置。计算所有粒子的目标函数，将其中的非劣解存在外部存档中，形成pareto解集；

step2 确定个体最优与群体最优

        由于是初始化过程，所以每个粒子的个体最优就是其本身；参考文献中群体最优的选取方式是采用密集距离进行排序，在初始化选取的时候为了方便，也可直接从pareto解集随机选取；

step3 更新惯性权重

        惯性权重是粒子群算法的一个重要参数，和算法的局部搜索能力以及全局搜索能力相关，在每次迭代时，将惯性权重按下面的公式更新：

 step4 更新粒子的速度和位置

         这一步和单目标粒子群算法一样：

 

 step5 交叉变异操作

        这一步也就是参考文献的“改进”所在，也是可以避免算法陷入局部最优的操作：

 step6 计算粒子的适应度并更新pareto最优解集

        参考文献中采用动态密集距离对pareto解集进行更新。简单来理解，就是判断一下pareto解集中各个解的距离，把距离较远的解留下(也就是保证解的分布不要太密集)，距离比较近的解淘汰，保证外部存档中pareto最优解的数量不超过上限。这一点也比较好理解，还是拿租房子来举例，假设现在pareto解集中有三个选项：A租金1800/月，面积100m²，B租金1600/月，面积80m²，C租金1799/月，面积99m²，三个解互相都是非支配的关系，但如果咱们外部存档最多只能存两个解，该淘汰谁呢？我们注意到，A和C之间两个目标的值都很接近，相似度很高，所以为了保证解的多样性，淘汰的时候就从A和C之中选一个，B可以保留下来(具体该淘汰谁没有具体算，只是举个例子方便理解)。这就是通过密集距离对pareto解集进行更新的原理，密集距离的计算公式如下：

step7 更新群体最优与个体最优

        群体最优和个体最优的更新方式也是MOPSO的研究重点，各种花里胡哨的方法都有，参考文献的方法就比较简单粗暴，随机选，具体方法如下：

 step8 迭代终止判断

        判断一下是否满足迭代次数要求，如果满足了就输出结果，不满足的话需要返回步骤3。

三、测试函数

        ZDT函数是多目标优化问题的常用的测试函数，具体公式如下：

1.ZDT1

​

2.ZDT2 

​

3.ZDT3 

​

4.ZDT4

​

四、部分代码

%% 清空环境
clc
clear
close all
%% 设置种群参数
sizepop = 200;                      % 初始种群个数
dim = 10;                           % 空间维数
ger = 100;                          % 最大迭代次数    
x_max=1*ones(1,dim);                % 位置上限
x_min=zeros(1,dim);                 % 位置下限
v_max=0.8*x_max;                    % 速度上限
v_min=-v_max;                       % 速度下限
w_start=0.9;                        % 惯性权重初始值
w_end=0.4;                          % 惯性权重结束值
c_1 = 1.5;                          % 自我学习因子
c_2 = 1.5;                          % 群体学习因子 
set_num_max=100;                    % pareto解集规模
%% 种群初始化
pop=x_min+rand(sizepop,dim).*(x_max-x_min);     % 初始化种群
pop_v=v_min+rand(sizepop,dim).*(v_max-v_min);   % 初始化种群速度                 
fitness=zeros(sizepop,2);                       % 所有个体的适应度
% 初始的适应度
for k=1:sizepop
    [a,b]=ZDT1(pop(k,:));
    fitness(k,:)=[a,b];
end
% 更新pareto解集
pareto_set=fitness(1,:);
pareto_pop=pop(1,:);
for k=2:sizepop
    nk=length(pareto_set(:,1));
    kk=1;
    index=0;% 判断是否放入pareto解集中的标志
    while kk<=nk
        % 如果pareto解集中被第k个个体支配就将其删除
        if fitness(k,:)<pareto_set(kk,:)
            pareto_set(kk,:)=[];
            pareto_pop(kk,:)=[];
            nk=nk-1;
            kk=kk-1;
        % 如果个体k被pareto解集中任意支配就退出循环
        elseif fitness(k,:)>pareto_set(kk,:)
            break;
        else
            index=index+1;
        end
        kk=kk+1;
    end
    if index>=length(pareto_set(:,1))
        pareto_set=[pareto_set;fitness(k,:)];
        pareto_pop=[pareto_pop;pop(k,:)];
    end
end
% 采用密集距离法更新全局最优和个体最优
I = crowd_dist(pareto_set);                     % 求所有个体的拥挤距离
% 采用逐一去除法更新pareto解集
while length(I)>set_num_max
    pareto_set(I==min(I),:)=[];
    I = crowd_dist(pareto_set);
end
% 在拥挤距离较大的前20%解中随机选择全局最优
if floor(0.2*length(I))<1
    r1=1;
else
    r1=randi([1,floor(0.2*length(I))]);
end
x_zbest=pareto_pop(r1,:);                       % 全局最优的位置
fitness_zbest=pareto_set(r1,:);                 % 种群的最优适应度
x_gbest=pop;                                    % 个体的最优位置   
fitness_gbest=fitness;                          % 个体的最优适应度
%% 迭代求pareto解集
iter=1;
while iter <= ger
    省略....
    iter=iter+1;
end

五、结果展示

        代码中以动画的形式更新pareto解集，这里只展示最终结果：

1.ZDT1

​

2.ZDT2

​

 3.ZDT3

 4.ZDT4

六、完整代码 

        完整代码获取方式参见链接：

多目标粒子群算法(MOPSO)的原理介绍以及matlab实现