NSGA-II原理及platEMO中Matlab代码实现

介绍

Non dominated sorting genetic algorithm -II

非支配排序遗传算法

①提出了快速非支配排序算法，一方面降低了计算的复杂度，另一方面它将父代种群跟子代种群进行合并，使得下一代的种群从双倍的空间中进行选取，从而保留了最为优秀的所有个体；
②引进精英策略，保证某些优良的种群个体在进化过程中不会被丢弃，从而提高了优化结果的精度；
③采用拥挤度和拥挤度比较算子，不但克服了NSGA中需要人为指定共享参数的缺陷，而且将其作为种群中个体间的比较标准，使得准Pareto域中的个体能均匀地扩展到整个Pareto域，保证了种群的多样性。

算法流程

伪代码如下：

NSGA-II的核心在于循环中：

①父代Pt通过选择、交叉、变异产生子代Qt，合并Pt和Qt为Rt
②对Rt进行非支配快速排序，按照等级排序为F1、F2……，为了保持种群的数量相等，依照等级高的进行选取，但会存在选取某一等级时，数量会超过种群的数量，为了保证种群的多样性，使用拥挤度计算方法来选取，拥挤度越高，选取的概率越大
③这样就产生了下一个父代Pt+1，再进入循环，直到达到迭代次数

选择快速排序得到的高等级和拥挤度大的个体的方法就是精英保留策略

非支配快速排序

首先要了解什么叫做支配
支配的定义如下：

$A(x1|y1)支配B(x2|y2)当：$
$$(x1<=x2andy1<=y2)and(x1<x2ory1<y2)$$

举一个例子：优化目标是同时最小化$f1$和$f2$，在上图中A和B的$f2$值是一样的，但A的$f1$值是小于B的$f1$值的，所以可以说A是支配B
matlab代码

all(x <= y) && any(x<y)

依据支配的定义，分别列出每个个体的被支配个数，以及支配的集合
对于被支配个数为0的个体（例如i），加入到Q（例如F1，表示等级rank为1）中，删除个体i，并对支配集合中有个体i的个体，将支配集合删除个体i，并支配个数-1。这类似于拓扑排序，以此往复，对每个个体进行等级排序。

拥挤度计算

将个体依据等级排序加入到新的父代Pt+1时，为了保持种群规模一致，需要对不能完全放完的某一等级的个体进行选择，选择的标准则是依据拥挤度，拥挤度越大，选择的概率越高。

个体的拥挤度是每个目标下的拥挤度总和
对于每个目标，先将种群按照目标值进行排序，对于最大和最小的目标值个体拥挤度设置为无穷大，其他的个体则依据相邻个体计算拥挤度

platEMO

platEMO是由安徽大学创建的一个平台，包含了很多优化算法，直接下载代码，放到matlab中即可

GUI界面更加直观

NSGA-II源码主要方法 Matlab实现

NSGA-II.m

classdef NSGAII < ALGORITHM
    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

EnvironmentalSelection.m

function [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Population,N)
    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;
    
    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);
    
    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;
    
    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
end

NDSort.m

function [FrontNo,MaxFNo] = NDSort(varargin)
    PopObj = varargin{1};
    [N,M]  = size(PopObj);%对传入的参数进行解析
    if nargin == 2
        nSort  = varargin{2};
    else 
        PopCon = varargin{2};%约束支配
        nSort  = varargin{3};
        Infeasible           = any(PopCon>0,2);
        PopObj(Infeasible,:) = repmat(max(PopObj,[],1),sum(Infeasible),1) + repmat(sum(max(0,PopCon(Infeasible,:)),2),1,M);
    end
    
    if M < 3 || N < 500
        % Use efficient non-dominated sort with sequential search (ENS-SS)
        [FrontNo,MaxFNo] = ENS_SS(PopObj,nSort);
    else
        %基于树形的排序优化（规模较大时）
        % Use tree-based efficient non-dominated sort (T-ENS)
        [FrontNo,MaxFNo] = T_ENS(PopObj,nSort);
    end
end

function [FrontNo,MaxFNo] = ENS_SS(PopObj,nSort)
    % unique类似于基数排序，找到不重复的,同时排/序
    [PopObj,~,Loc] = unique(PopObj,'rows');
    Table   = hist(Loc,1:max(Loc));
    [N,M]   = size(PopObj);
    FrontNo = inf(1,N);
    MaxFNo  = 0;
    % while循环控制层数
    while sum(Table(FrontNo<inf)) < min(nSort,length(Loc))
        MaxFNo = MaxFNo + 1;
        % 第一层循环，遍历所有个体，找到是否有属于同一个非支配层的解
        for i = 1 : N
            if FrontNo(i) == inf
                Dominated = false;
                % 第二层循环，分析j和已经访问过的解是否存在支配关系，如果存在支配关系
                % break；否则，说明j属性MaxFNo对应的层
                for j = i-1 : -1 : 1
                    if FrontNo(j) == MaxFNo
                        m = 2;
                        % 比较所有目标，如果当前被比较的个体i，被个体j所支配，m自增1
                        while m <= M && PopObj(i,m) >= PopObj(j,m)
                            m = m + 1;
                        end
                        % 被支配 设置Dominated 为true
                        Dominated = m > M;
                        if Dominated || M == 2
                            break;
                        end
                    end
                end
                % 如果没有任何的个体支配当前个体i，则设置MaxFNo
                if ~Dominated
                    FrontNo(i) = MaxFNo;
                end
            end
        end
    end
    % 放回到原始数据中
    FrontNo = FrontNo(:,Loc);
end

CrowdingDistance.m

function CrowdDis = CrowdingDistance(PopObj,FrontNo)
    [N,M] = size(PopObj);
    if nargin < 2
        FrontNo = ones(1,N);
    end
    CrowdDis = zeros(1,N);
    Fronts   = setdiff(unique(FrontNo),inf);
    for f = 1 : length(Fronts)
        Front = find(FrontNo==Fronts(f));
        Fmax  = max(PopObj(Front,:),[],1);
        Fmin  = min(PopObj(Front,:),[],1);
        for i = 1 : M
            [~,Rank] = sortrows(PopObj(Front,i));
            CrowdDis(Front(Rank(1)))   = inf;
            CrowdDis(Front(Rank(end))) = inf;
            for j = 2 : length(Front)-1
                CrowdDis(Front(Rank(j))) = CrowdDis(Front(Rank(j)))+(PopObj(Front(Rank(j+1)),i)-PopObj(Front(Rank(j-1)),i))/(Fmax(i)-Fmin(i));
            end
        end
    end
end

参考

【加中文字幕】理解NSGA-Ⅱ的工作原理】
【PlatEMO指北：P7-NSGA-II完全解析（下）源码分析】