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

一般遗传算法的流程：

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

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

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

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

A、快速非支配排序方法

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

算法流程如下：

首先计算每个解的np和Sp，这需要O(MN2)时间复杂度，同时得到了第一前沿面

F

1

\mathcal{F1}

F1。第二部分就是计算其余的前沿面，需要遍历前沿面中的每个解以及这个解的Sp集合，将对应解p的np值减一，如果np值减为0了，就加入下一前沿面集合，这部分需要O(N2)时间复杂度。由于第一部分和第二部分是分开的，所以总的时间复杂度是O(MN2)。

B、拥挤比较算子

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

具体算法流程如下：

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

C、主程序

• 种群初始化：随机创建父代种群P0 ，为父代种群利用快速非支配排序算法计算出非支配等级，并利用二进制锦标赛机制重组以及变异算子生成大小为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(i2-1,:)=(population(k(1),:)+population(k(2),:))/2+beta.(population(k(1),:)-population(k(2),:))./2; %产生第一个子代
  newpopulation(i2,:)=(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')

结果