多目标优化经典算法——NSGA-II

参考博客链接
https://blog.csdn.net/qq_35414569/article/details/79639848?utm_medium=distribute.pc_relevant.none-task-blog-2_defaultbaidujs_baidulandingword~default-0.no_search_link&spm=1001.2101.3001.4242
因为NSGA-II算法是一种遗传算法，所以首先搞清楚遗传算法的流程。

遗传算法流程

一般遗传算法的流程：

1. 种群初始化
2. 计算每个个体的适应度
3. 选择
4. 交叉
5. 变异

根据是否满足解的精度要求和迭代次数来判断是否进行下一轮的遗传进化。

NSGA算法存在的3个问题

1. O(MN^3)计算时间复杂度（其中M代表目标个数，N代表种群个数）
2. 非精英机制方法
3. 需要指定一个共享参数

NSGA-II算法

NSGA-II算法主要由以下三个部分组成
A、快速非支配排序方法
B、拥挤比较算子
C、主程序

A、快速非支配排序方法

• 传统排序方法：时间复杂度O(MN³)，M是目标个数，N是种群个数。为了计算第一非支配前沿面，需要判断每个解和种群中的其他解的支配关系。一个解和其他解的支配关系需要O(MN)复杂度，每个解和其他解的支配关系需要O(MN²)复杂度。最坏的情况，N个解各自构成一个前沿面，时间复杂度就是O(MN³)。
• 快速非支配排序方法：首先需要计算两个变量 1）支配计数n_p，即支配解p的解的数量。2）S_p，即解p支配的解的集合。
  
  如上图所示，对解c而言，解c被a,b两个解支配，所以n_c为2。对解b而言，解b支配c,d,e，所以S_b={c,d,e}。

算法流程如下：

首先计算每个解的n_p和S_p，这需要O(MN²)时间复杂度，同时得到了第一前沿面$\mathcal{F}1$。第二部分就是计算其余的前沿面，需要遍历前沿面中的每个解以及这个解的S_p集合，将对应解p的n_p值减一，如果n_p值减为0了，就加入下一前沿面集合，这部分需要O(N²)时间复杂度。由于第一部分和第二部分是分开的，所以总的时间复杂度是O(MN²)。

B、拥挤比较算子

密度估计：对同一前沿面的解按照目标函数值排序，计算每个解在该目标函数下两侧解的目标函数归一化差值，然后将所有目标函数下的分量累加作为拥挤系数。

具体算法流程如下：

如图所示，直观上看，解i的拥挤系数就是虚线矩形周长的一半。

C、主程序

• 种群初始化：随机创建父代种群P₀ ，为父代种群利用快速非支配排序算法计算出非支配等级，并利用二进制锦标赛机制重组以及变异算子生成大小为N的子代种群。
• 子代个体选择
  
  代码实现
  参考博客链接
  https://www.omegaxyz.com/2017/05/04/nsga2matlabzdt1/
  主体代码是这篇博客的，仅将非支配排序部分修改为了原始论文中的

%% NSGA-II算法
clear; close all; clc; tic;

%% 初始化/参数设定
generations=100;                                %迭代次数
popnum=100;                                     %种群大小(须为偶数)
poplength=30;                                   %个体长度
minvalue=repmat(zeros(1,poplength),popnum,1);   %个体最小值
maxvalue=repmat(ones(1,poplength),popnum,1);    %个体最大值    
%% 随机创建父代种群  
population=rand(popnum,poplength).*(maxvalue-minvalue)+minvalue;    %产生新的初始种群

%% 开始迭代进化
for gene=1:generations                      %开始迭代
    %% 交叉
    newpopulation=zeros(popnum,poplength);  %子代种群
    for i=1:popnum/2                        %交叉产生子代
            k=randperm(popnum);                 %从种群中随机选出两个父母,不采用二进制联赛方法
            beta=(-1).^round(rand(1,poplength)).*abs(randn(1,poplength))*1.481; %采用正态分布交叉产生两个子代
            newpopulation(i*2-1,:)=(population(k(1),:)+population(k(2),:))/2+beta.*(population(k(1),:)-population(k(2),:))./2;    %产生第一个子代           
            newpopulation(i*2,:)=(population(k(1),:)+population(k(2),:))/2-beta.*(population(k(1),:)-population(k(2),:))./2;      %产生第二个子代
    end
    
    %% 变异
    k=rand(size(newpopulation));    %随机选定要变异的基因位
    miu=rand(size(newpopulation));  %采用多项式变异
    temp=k<1/poplength & miu<0.5;   %要变异的基因位
    newpopulation(temp)=newpopulation(temp)+(maxvalue(temp)-minvalue(temp)).*((2.*miu(temp)+(1-2.*miu(temp)).*(1-(newpopulation(temp)-minvalue(temp))./(maxvalue(temp)-minvalue(temp))).^21).^(1/21)-1);        %变异情况一
    newpopulation(temp)=newpopulation(temp)+(maxvalue(temp)-minvalue(temp)).*(1-(2.*(1-miu(temp))+2.*(miu(temp)-0.5).*(1-(maxvalue(temp)-newpopulation(temp))./(maxvalue(temp)-minvalue(temp))).^21).^(1/21));  %变异情况二
    
    %% 越界处理/种群合并
    newpopulation(newpopulation>maxvalue)=maxvalue(newpopulation>maxvalue); %子代越上界处理
    newpopulation(newpopulation<minvalue)=minvalue(newpopulation<minvalue); %子代越下界处理
    newpopulation=[population;newpopulation];                               %合并父子种群
    
    %% 计算目标函数值
    functionvalue=zeros(size(newpopulation,1),2);           %合并后种群的各目标函数值,这里的问题是ZDT1
    functionvalue(:,1)=newpopulation(:,1);                  %计算第一维目标函数值
    g=1+9*sum(newpopulation(:,2:poplength),2)./(poplength-1);
    functionvalue(:,2)=g.*(1-(newpopulation(:,1)./g).^0.5); %计算第二维目标函数值
    
    %% 非支配排序，NSGA-II论文中的算法
    Sp = zeros(size(newpopulation,1));              %Sp解p支配的解的集合
    np = zeros(size(newpopulation,1),1);            %支配解p的解的数量
    Prank = zeros(size(newpopulation,1),1);         %每个解的前沿面编号
    F = zeros(size(newpopulation,1));               %保存每层前沿面的个体
    for p = 1:size(newpopulation,1)
        for q = 1:size(newpopulation,1)
            if isDomin(functionvalue,p,q)
                Sp(p,sum(Sp(p,:)~=0)+1) = q;
            elseif isDomin(functionvalue,q,p)
                np(p) = np(p)+1; 
            end
        end
        if np(p)==0
            Prank(p) = 1;
            F(1,sum(F(1,:)~=0)+1) = p;
        end
    end
    i = 1;
    while sum(F(i,:)~=0)
        Q = zeros(1,size(F,2));
        for p = 1:sum(F(i,:)~=0)
            for q = 1:sum(Sp(F(i,p),:)~=0)
                np(Sp(F(i,p),q)) = np(Sp(F(i,p),q)) - 1;
                if np(Sp(F(i,p),q))==0
                    Prank(Sp(F(i,p),q)) = i+1;
                    Q(sum(Q(:)~=0)+1) = Sp(F(i,p),q);
                end
            end
        end
        i = i+1;
        F(i,:) = Q;
    end
    
    %% 计算拥挤距离/选出下一代个体  
    fnum = 0;
    psum = 0;
    while psum <= popnum                        %判断前多少个面的个体能完全放入下一代种群
        fnum = fnum + 1;
        psum = psum + sum(F(fnum,:)~=0);
    end
    num = popnum - (psum - sum(F(fnum,:)~=0));  %num表示在当前前沿面需要多少个个体
    population(1:popnum-num,:) = newpopulation(Prank<fnum,:);
    % 计算当前前沿面个体的拥挤距离
    I = F(fnum,:);
    l = sum(I(:)~=0);
    I = I(1:l);
    I_distance = zeros(l,1);
    fmax = max(functionvalue(I,:),[],1);
    fmin = min(functionvalue(I,:),[],1);
    for m = 1:size(functionvalue,2)
       [~,oldsite] = sortrows(functionvalue(I,m));
       I_distance(oldsite(1)) = inf;
       I_distance(oldsite(end)) = inf;
       for i = 2:l-1
           I_distance(oldsite(i)) = I_distance(oldsite(i)) + (functionvalue(I(oldsite(i+1)),m)-functionvalue(I(oldsite(i-1)),m))/(fmax(m)-fmin(m));
       end
    end
    I = sortrows([I_distance,I'],1,'descend');                          %按拥挤距离降序排序
    population(popnum-num+1:popnum,:) = newpopulation(I(1:num,2),:);    %将第fnum个面上拥挤距离较大的前num个个体复制入下一代 
end
    
%% 程序输出
fprintf('已完成,耗时%4s秒\n',num2str(toc));          %程序最终耗时
output=sortrows(functionvalue(Prank==1,:));    %最终结果:种群中非支配解的函数值
plot(output(:,1),output(:,2),'*b');                 %作图
axis([0,1,0,1]);xlabel('F_1');ylabel('F_2');title('ZDT1')    
    

结果