zoukankan      html  css  js  c++  java
  • NSGA3理解(NSGA3算法及其MATLAB版本实现)

    NSGA3算法及其MATLAB版本实现

    NSGA3理解

    对非支配遗传算法三的一些理解

    GitHub上的C++版本:https://github.com/adanjoga/nsga3

     

    晓风从platEMO中摘出来的NSGA-III算法代码:

    NSGAIII_main.m

     

     1 clc,clear
     2 global N M D  PopCon name gen
     3 
     4 N = 400;                        % 种群个数
     5 M = 3;                          % 目标个数
     6 name = 'DTLZ1';                 % 测试函数选择,目前只有:DTLZ1、DTLZ2、DTLZ3
     7 gen = 500;                      %迭代次数
     8 %% Generate the reference points and random population
     9 [Z,N] = UniformPoint(N,M);        % 生成一致性参考解
    10 [res,Population,PF] = funfun(); % 生成初始种群与目标值
    11 Pop_objs = CalObj(Population); % 计算适应度函数值
    12 Zmin  = min(Pop_objs(all(PopCon<=0,2),:),[],1); %求理想点,其实PopCon是处理有约束问题的,这里并没有用到
    13 
    14 %% Optimization
    15 for i = 1:gen
    16     MatingPool = TournamentSelection(2,N,sum(max(0,PopCon),2));
    17     Offspring  = GA(Population(MatingPool,:));
    18     Offspring_objs = CalObj(Offspring);
    19     Zmin       = min([Zmin;Offspring_objs],[],1);
    20     Population = EnvironmentalSelection([Population;Offspring],N,Z,Zmin);
    21     Popobj = CalObj(Population);
    22     if(M<=3)
    23         plot3(Popobj(:,1),Popobj(:,2),Popobj(:,3),'ro')
    24         title(num2str(i));
    25         drawnow
    26     end
    27 end
    28 if(M<=3)
    29     hold on
    30     plot3(PF(:,1),PF(:,2),PF(:,3),'g*')
    31 else
    32     for i=1:size(Popobj,1)
    33         plot(Popobj(i,:))
    34         hold on
    35     end
    36 end
    37 %%IGD
    38 score = IGD(Popobj,PF)

     

     

    UniformPoint.m

     1 function [W,N] = UniformPoint(N,M)
     2 %UniformPoint - Generate a set of uniformly distributed points on the unit
     3 %hyperplane.
     4 %
     5 %   [W,N] = UniformPoint(N,M) returns approximately N uniformly distributed
     6 %   points with M objectives on the unit hyperplane.
     7 %
     8 %   Due to the requirement of uniform distribution, the number of points
     9 %   cannot be arbitrary, and the number of points in W may be slightly
    10 %   smaller than the predefined size N.
    11 %
    12 %   Example:
    13 %       [W,N] = UniformPoint(275,10)
    14     H1 = 1;
    15     while nchoosek(H1+M-1,M-1) <= N
    16         H1 = H1 + 1;
    17     end
    18     H1=H1 - 1;
    19     W = nchoosek(1:H1+M-1,M-1) - repmat(0:M-2,nchoosek(H1+M-1,M-1),1) - 1;%nchoosek(v,M),v是一个向量,返回一个矩阵,该矩阵的行由每次从v中的M个元素取出k个取值的所有可能组合构成。比如:v=[1,2,3],n=2,输出[1,2;1,3;2,3]
    20     W = ([W,zeros(size(W,1),1)+H1]-[zeros(size(W,1),1),W])/H1;
    21     if H1 < M
    22         H2 = 0;
    23         while nchoosek(H1+M-1,M-1)+nchoosek(H2+M-1,M-1) <= N
    24             H2 = H2 + 1;
    25         end
    26         H2 = H2 - 1;
    27         if H2 > 0
    28             W2 = nchoosek(1:H2+M-1,M-1) - repmat(0:M-2,nchoosek(H2+M-1,M-1),1) - 1;
    29             W2 = ([W2,zeros(size(W2,1),1)+H2]-[zeros(size(W2,1),1),W2])/H2;
    30             W  = [W;W2/2+1/(2*M)];
    31         end
    32     end
    33     W = max(W,1e-6);
    34     N = size(W,1);
    35 end

    funfun.m

     1 function [PopObj,PopDec,P] =  funfun()
     2 
     3 global M D lower upper encoding N PopCon cons cons_flag name
     4 
     5     %% Initialization
     6     D = M + 4;
     7     lower    = zeros(1,D);
     8     upper    = ones(1,D);
     9     encoding = 'real';
    10     switch encoding
    11         case 'binary'
    12             PopDec = randi([0,1],N,D);
    13         case 'permutation'
    14             [~,PopDec] = sort(rand(N,D),2);
    15         otherwise
    16             PopDec = unifrnd(repmat(lower,N,1),repmat(upper,N,1));
    17     end
    18     cons = zeros(size(PopDec,1),1);
    19     cons_flag = 1;
    20     PopCon = cons;
    21     %% Calculate objective values
    22     switch name
    23         case 'DTLZ1'
    24             g      = 100*(D-M+1+sum((PopDec(:,M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,M:end)-0.5)),2));
    25             PopObj = 0.5*repmat(1+g,1,M).*fliplr(cumprod([ones(N,1),PopDec(:,1:M-1)],2)).*[ones(N,1),1-PopDec(:,M-1:-1:1)];
    26             P = UniformPoint(N,M)/2;
    27         case 'DTLZ2'
    28             g      = sum((PopDec(:,M:end)-0.5).^2,2);
    29             PopObj = repmat(1+g,1,M).*fliplr(cumprod([ones(size(g,1),1),cos(PopDec(:,1:M-1)*pi/2)],2)).*[ones(size(g,1),1),sin(PopDec(:,M-1:-1:1)*pi/2)];
    30             P = UniformPoint(N,M);
    31             P = P./repmat(sqrt(sum(P.^2,2)),1,M);
    32         case 'DTLZ3'
    33             g      = 100*(D-M+1+sum((PopDec(:,M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,M:end)-0.5)),2));
    34             PopObj = repmat(1+g,1,M).*fliplr(cumprod([ones(N,1),cos(PopDec(:,1:M-1)*pi/2)],2)).*[ones(N,1),sin(PopDec(:,M-1:-1:1)*pi/2)];
    35             P = UniformPoint(N,M);
    36             P = P./repmat(sqrt(sum(P.^2,2)),1,M);
    37     end
    38    %% Sample reference points on Pareto front
    39     %P = UniformPoint(N,obj.Global.M)/2;
    40 end

    CalObj.m

     1 function PopObj = CalObj(PopDec)
     2 global M D name
     3     NN = size(PopDec,1);
     4     switch name
     5         case 'DTLZ1'
     6             g      = 100*(D-M+1+sum((PopDec(:,M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,M:end)-0.5)),2));
     7             PopObj = 0.5*repmat(1+g,1,M).*fliplr(cumprod([ones(NN,1),PopDec(:,1:M-1)],2)).*[ones(NN,1),1-PopDec(:,M-1:-1:1)];
     8         case 'DTLZ2'
     9             g      = sum((PopDec(:,M:end)-0.5).^2,2);
    10             PopObj = repmat(1+g,1,M).*fliplr(cumprod([ones(size(g,1),1),cos(PopDec(:,1:M-1)*pi/2)],2)).*[ones(size(g,1),1),sin(PopDec(:,M-1:-1:1)*pi/2)];
    11         case 'DTLZ3'
    12             g      = 100*(D-M+1+sum((PopDec(:,M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,M:end)-0.5)),2));
    13             PopObj = repmat(1+g,1,M).*fliplr(cumprod([ones(NN,1),cos(PopDec(:,1:M-1)*pi/2)],2)).*[ones(NN,1),sin(PopDec(:,M-1:-1:1)*pi/2)];
    14     end
    15 end

    TournamentSelection.m

     1 function index = TournamentSelection(K,N,varargin)
     2 %TournamentSelection - Tournament selection.
     3 %
     4 %   P = TournamentSelection(K,N,fitness1,fitness2,...) returns the indices
     5 %   of N solutions by K-tournament selection based on their fitness values.
     6 %   In each selection, the candidate having the MINIMUM fitness1 value will
     7 %   be selected; if more than one candidates have the same minimum value of
     8 %   fitness1, then compare their fitness2 values, and so on.
     9 %
    10 %   Example:
    11 %       P = TournamentSelection(2,100,FrontNo)
    12 
    13 %------------------------------- Copyright --------------------------------
    14 % Copyright (c) 2018-2019 BIMK Group. You are free to use the PlatEMO for
    15 % research purposes. All publications which use this platform or any code
    16 % in the platform should acknowledge the use of "PlatEMO" and reference "Ye
    17 % Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
    18 % for evolutionary multi-objective optimization [educational forum], IEEE
    19 % Computational Intelligence Magazine, 2017, 12(4): 73-87".
    20 %--------------------------------------------------------------------------
    21 
    22     varargin = cellfun(@(S) reshape(S,[],1),varargin,'UniformOutput',false);
    23     [~,rank] = sortrows([varargin{:}]);
    24     [~,rank] = sort(rank);
    25     Parents  = randi(length(varargin{1}),K,N);
    26     [~,best] = min(rank(Parents),[],1);
    27     index    = Parents(best+(0:N-1)*K);
    28 end

    GA.m

      1 function Offspring = GA(Parent,Parameter)
      2 global encoding lower upper
      3 %GA - Genetic operators for real, binary, and permutation based encodings.
      4 %
      5 %   Off = GA(P) returns the offsprings generated by genetic operators,
      6 %   where P1 is a set of parents. If P is an array of INDIVIDUAL objects,
      7 %   then Off is also an array of INDIVIDUAL objects; while if P is a matrix
      8 %   of decision variables, then Off is also a matrix of decision variables,
      9 %   i.e., the offsprings will not be evaluated. P is split into two subsets
     10 %   P1 and P2 with the same size, and each object/row of P1 and P2 is used
     11 %   to generate TWO offsprings. Different operators are used for real,
     12 %   binary, and permutation based encodings, respectively.
     13 %
     14 %   Off = GA(P,{proC,disC,proM,disM}) specifies the parameters of
     15 %   operators, where proC is the probabilities of doing crossover, disC is
     16 %   the distribution index of simulated binary crossover, proM is the
     17 %   expectation of number of bits doing mutation, and disM is the
     18 %   distribution index of polynomial mutation.
     19 %
     20 %   Example:
     21 %       Off = GA(Parent)
     22 %       Off = GA(Parent.decs,{1,20,1,20})
     23 %
     24 %   See also GAhalf
     25 
     26 %------------------------------- Reference --------------------------------
     27 % [1] K. Deb, K. Sindhya, and T. Okabe, Self-adaptive simulated binary
     28 % crossover for real-parameter optimization, Proceedings of the 9th Annual
     29 % Conference on Genetic and Evolutionary Computation, 2007, 1187-1194.
     30 % [2] K. Deb and M. Goyal, A combined genetic adaptive search (GeneAS) for
     31 % engineering design, Computer Science and informatics, 1996, 26: 30-45.
     32 % [3] L. Davis, Applying adaptive algorithms to epistatic domains,
     33 % Proceedings of the International Joint Conference on Artificial
     34 % Intelligence, 1985, 162-164.
     35 % [4] D. B. Fogel, An evolutionary approach to the traveling salesman
     36 % problem, Biological Cybernetics, 1988, 60(2): 139-144.
     37 %------------------------------- Copyright --------------------------------
     38 % Copyright (c) 2018-2019 BIMK Group. You are free to use the PlatEMO for
     39 % research purposes. All publications which use this platform or any code
     40 % in the platform should acknowledge the use of "PlatEMO" and reference "Ye
     41 % Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
     42 % for evolutionary multi-objective optimization [educational forum], IEEE
     43 % Computational Intelligence Magazine, 2017, 12(4): 73-87".
     44 %--------------------------------------------------------------------------
     45 
     46     %% Parameter setting
     47     if nargin > 1
     48         [proC,disC,proM,disM] = deal(Parameter{:});
     49     else
     50         [proC,disC,proM,disM] = deal(1,20,1,20);
     51     end
     52     if isa(Parent(1),'INDIVIDUAL')
     53         calObj = true;
     54         Parent = Parent.decs;
     55     else
     56         calObj = false;
     57     end
     58     Parent1 = Parent(1:floor(end/2),:);
     59     Parent2 = Parent(floor(end/2)+1:floor(end/2)*2,:);
     60     [N,D]   = size(Parent1);
     61  
     62     switch encoding
     63         case 'binary'
     64             %% Genetic operators for binary encoding
     65             % One point crossover
     66             k = repmat(1:D,N,1) > repmat(randi(D,N,1),1,D);
     67             k(repmat(rand(N,1)>proC,1,D)) = false;
     68             Offspring1    = Parent1;
     69             Offspring2    = Parent2;
     70             Offspring1(k) = Parent2(k);
     71             Offspring2(k) = Parent1(k);
     72             Offspring     = [Offspring1;Offspring2];
     73             % Bitwise mutation
     74             Site = rand(2*N,D) < proM/D;
     75             Offspring(Site) = ~Offspring(Site);
     76         case 'permutation'
     77             %% Genetic operators for permutation based encoding
     78             % Order crossover
     79             Offspring = [Parent1;Parent2];
     80             k = randi(D,1,2*N);
     81             for i = 1 : N
     82                 Offspring(i,k(i)+1:end)   = setdiff(Parent2(i,:),Parent1(i,1:k(i)),'stable');
     83                 Offspring(i+N,k(i)+1:end) = setdiff(Parent1(i,:),Parent2(i,1:k(i)),'stable');
     84             end
     85             % Slight mutation
     86             k = randi(D,1,2*N);
     87             s = randi(D,1,2*N);
     88             for i = 1 : 2*N
     89                 if s(i) < k(i)
     90                     Offspring(i,:) = Offspring(i,[1:s(i)-1,k(i),s(i):k(i)-1,k(i)+1:end]);
     91                 elseif s(i) > k(i)
     92                     Offspring(i,:) = Offspring(i,[1:k(i)-1,k(i)+1:s(i)-1,k(i),s(i):end]);
     93                 end
     94             end
     95         otherwise
     96             %% Genetic operators for real encoding
     97             % Simulated binary crossover
     98             beta = zeros(N,D);
     99             mu   = rand(N,D);
    100             beta(mu<=0.5) = (2*mu(mu<=0.5)).^(1/(disC+1));
    101             beta(mu>0.5)  = (2-2*mu(mu>0.5)).^(-1/(disC+1));
    102             beta = beta.*(-1).^randi([0,1],N,D);
    103             beta(rand(N,D)<0.5) = 1;
    104             beta(repmat(rand(N,1)>proC,1,D)) = 1;
    105             Offspring = [(Parent1+Parent2)/2+beta.*(Parent1-Parent2)/2
    106                          (Parent1+Parent2)/2-beta.*(Parent1-Parent2)/2];
    107             % Polynomial mutation
    108             Lower = repmat(lower,2*N,1);
    109             Upper = repmat(upper,2*N,1);
    110             Site  = rand(2*N,D) < proM/D;
    111             mu    = rand(2*N,D);
    112             temp  = Site & mu<=0.5;
    113             Offspring       = min(max(Offspring,Lower),Upper);
    114             Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
    115                               (1-(Offspring(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    116             temp = Site & mu>0.5; 
    117             Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
    118                               (1-(Upper(temp)-Offspring(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
    119     end
    120     if calObj
    121         Offspring = INDIVIDUAL(Offspring);
    122     end
    123 end

    EnvironmentalSelection.m

     1 function Population = EnvironmentalSelection(Population,N,Z,Zmin)
     2 global PopCon
     3 % The environmental selection of NSGA-III
     4 
     5 %------------------------------- Copyright --------------------------------
     6 % Copyright (c) 2018-2019 BIMK Group. You are free to use the PlatEMO for
     7 % research purposes. All publications which use this platform or any code
     8 % in the platform should acknowledge the use of "PlatEMO" and reference "Ye
     9 % Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
    10 % for evolutionary multi-objective optimization [educational forum], IEEE
    11 % Computational Intelligence Magazine, 2017, 12(4): 73-87".
    12 %--------------------------------------------------------------------------
    13 
    14     if isempty(Zmin)
    15         Zmin = ones(1,size(Z,2));
    16     end
    17 
    18     %% Non-dominated sorting
    19     Population_objs = CalObj(Population);
    20     [FrontNo,MaxFNo] = NDSort(Population_objs,PopCon,N);
    21     Next = FrontNo < MaxFNo;
    22     
    23     %% Select the solutions in the last front
    24     Last   = find(FrontNo==MaxFNo);
    25     Choose = LastSelection(Population_objs(Next,:),Population_objs(Last,:),N-sum(Next),Z,Zmin);
    26     Next(Last(Choose)) = true;
    27     % Population for next generation
    28     Population = Population(Next,:);
    29 end
    30 
    31 function Choose = LastSelection(PopObj1,PopObj2,K,Z,Zmin)
    32 % Select part of the solutions in the last front
    33 
    34     PopObj = [PopObj1;PopObj2] - repmat(Zmin,size(PopObj1,1)+size(PopObj2,1),1);
    35     [N,M]  = size(PopObj);
    36     N1     = size(PopObj1,1);
    37     N2     = size(PopObj2,1);
    38     NZ     = size(Z,1);
    39 
    40     %% Normalization
    41     % Detect the extreme points
    42     Extreme = zeros(1,M);
    43     w       = zeros(M)+1e-6+eye(M);
    44     for i = 1 : M
    45         [~,Extreme(i)] = min(max(PopObj./repmat(w(i,:),N,1),[],2));
    46     end
    47     % Calculate the intercepts of the hyperplane constructed by the extreme
    48     % points and the axes
    49     Hyperplane = PopObj(Extreme,:)ones(M,1);%A的逆矩阵是 A1,则 B/A = B*A1,AB = A1*B
    50     a = 1./Hyperplane;
    51     if any(isnan(a))%若a的元素为NaN(非数值),在对应位置上返回逻辑1(真),否则返回逻辑0(假)
    52         a = max(PopObj,[],1)';
    53     end
    54     % Normalization
    55     %a = a-Zmin';
    56     PopObj = PopObj./repmat(a',N,1);%如果a、b是矩阵,a./b就是a、b中对应的每个元素相除,得到一个新的矩阵;
    57     
    58     %% Associate each solution with one reference point
    59     % Calculate the distance of each solution to each reference vector
    60     Cosine   = 1 - pdist2(PopObj,Z,'cosine');%cosine’ 针对向量,夹角余弦距离Cosine distance(‘cosine’)余弦相似度= 1-余弦距离 
    61     
    62     %杰卡德距离Jaccard distance(‘jaccard’)Jaccard距离常用来处理仅包含非对称的二元(0-1)属性的对象。很显然,Jaccard距离不关心0-0匹配[1]。
    63     %夹角余弦距离Cosine distance(‘cosine’)与Jaccard距离相比,Cosine距离不仅忽略0-0匹配,而且能够处理非二元向量,即考虑到变量值的大小。
    64     %对这两者,距离与相似度和为一。
    65     %https://www.cnblogs.com/chaosimple/archive/2013/06/28/3160839.html
    66 
    67     Distance = repmat(sqrt(sum(PopObj.^2,2)),1,NZ).*sqrt(1-Cosine.^2);
    68     %Distance = pdist2(PopObj,Z,'euclidean');%上面注释的两个语句也可以哦
    69     % Associate each solution with its nearest reference point
    70     [d,pi] = min(Distance',[],1);
    71 
    72     %% Calculate the number of associated solutions except for the last front of each reference point
    73     rho = hist(pi(1:N1),1:NZ);%除了最后front的每一个参考点,计算关联解的个数
    74     
    75     %% Environmental selection
    76     Choose  = false(1,N2);
    77     Zchoose = true(1,NZ);
    78     % Select K solutions one by one
    79     while sum(Choose) < K
    80         % Select the least crowded reference point
    81         Temp = find(Zchoose);
    82         Jmin = find(rho(Temp)==min(rho(Temp)));%所有参考点中与之关联个数最少的那个参考点
    83         j    = Temp(Jmin(randi(length(Jmin))));
    84         I    = find(Choose==0 & pi(N1+1:end)==j);
    85         % Then select one solution associated with this reference point
    86         if ~isempty(I)
    87             if rho(j) == 0
    88                 [~,s] = min(d(N1+I));
    89             else
    90                 s = randi(length(I));
    91             end
    92             Choose(I(s)) = true;
    93             rho(j) = rho(j) + 1;
    94         else
    95             Zchoose(j) = false;
    96         end
    97     end
    98 end

     NDSort.m

      1 function [FrontNo,MaxFNo] = NDSort(varargin)
      2 %NDSort - Do non-dominated sorting by efficient non-dominated sort.
      3 %
      4 %   FrontNo = NDSort(F,s) does non-dominated sorting on F, where F is the
      5 %   matrix of objective values of a set of individuals, and s is the number
      6 %   of individuals to be sorted at least. FrontNo(i) denotes the front
      7 %   number of the i-th individual. The individuals have not been sorted are
      8 %   assigned a front number of inf.
      9 %
     10 %   FrontNo = NDSort(F,C,s) does non-dominated sorting based on constrained
     11 %   domination, where C is the matrix of constraint values of the
     12 %   individuals. In this case, feasible solutions always dominate
     13 %   infeasible solutions, and one infeasible solution dominates another
     14 %   infeasible solution if the former has a smaller overall constraint
     15 %   violation than the latter.
     16 %
     17 %   In particular, s = 1 indicates finding only the first non-dominated
     18 %   front, s = size(F,1)/2 indicates sorting only half the population
     19 %   (which is often used in the algorithm), and s = inf indicates sorting
     20 %   the whole population.
     21 %
     22 %   [FrontNo,K] = NDSort(...) also returns the maximum front number besides
     23 %   inf.
     24 %
     25 %   Example:
     26 %       [FrontNo,MaxFNo] = NDSort(PopObj,1)
     27 %       [FrontNo,MaxFNo] = NDSort(PopObj,PopCon,inf)
     28 
     29 %------------------------------- Reference --------------------------------
     30 % [1] X. Zhang, Y. Tian, R. Cheng, and Y. Jin, An efficient approach to
     31 % nondominated sorting for evolutionary multiobjective optimization, IEEE
     32 % Transactions on Evolutionary Computation, 2015, 19(2): 201-213.
     33 % [2] X. Zhang, Y. Tian, R. Cheng, and Y. Jin, A decision variable
     34 % clustering based evolutionary algorithm for large-scale many-objective
     35 % optimization, IEEE Transactions on Evolutionary Computation, 2018, 22(1):
     36 % 97-112.
     37 %------------------------------- Copyright --------------------------------
     38 % Copyright (c) 2018-2019 BIMK Group. You are free to use the PlatEMO for
     39 % research purposes. All publications which use this platform or any code
     40 % in the platform should acknowledge the use of "PlatEMO" and reference "Ye
     41 % Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
     42 % for evolutionary multi-objective optimization [educational forum], IEEE
     43 % Computational Intelligence Magazine, 2017, 12(4): 73-87".
     44 %--------------------------------------------------------------------------
     45 
     46     PopObj = varargin{1};
     47     [N,M]  = size(PopObj);
     48     if nargin == 2
     49         nSort  = varargin{2};
     50     else
     51         PopCon = varargin{2};
     52         nSort  = varargin{3};
     53         Infeasible           = any(PopCon>0,2);
     54         PopObj(Infeasible,:) = repmat(max(PopObj,[],1),sum(Infeasible),1) + repmat(sum(max(0,PopCon(Infeasible,:)),2),1,M);
     55     end
     56     if M < 5 || N < 500
     57         % Use efficient non-dominated sort with sequential search (ENS-SS)
     58         [FrontNo,MaxFNo] = ENS_SS(PopObj,nSort);
     59     else
     60         % Use tree-based efficient non-dominated sort (T-ENS)
     61         [FrontNo,MaxFNo] = T_ENS(PopObj,nSort);
     62     end
     63 end
     64 
     65 function [FrontNo,MaxFNo] = ENS_SS(PopObj,nSort)
     66     [PopObj,~,Loc] = unique(PopObj,'rows');
     67     Table          = hist(Loc,1:max(Loc));
     68     [N,M]          = size(PopObj);
     69     [PopObj,rank]  = sortrows(PopObj);
     70 %     if ~(all(PopObj(:) == PopObj1(:)))
     71 %         disp(2)
     72 %     end
     73     FrontNo        = inf(1,N);
     74     MaxFNo         = 0;
     75     while sum(Table(FrontNo<inf)) < min(nSort,length(Loc))
     76         MaxFNo = MaxFNo + 1;
     77         for i = 1 : N
     78             if FrontNo(i) == inf
     79                 Dominated = false;
     80                 for j = i-1 : -1 : 1
     81                     if FrontNo(j) == MaxFNo
     82                         m = 2;
     83                         while m <= M && PopObj(i,m) >= PopObj(j,m)
     84                             m = m + 1;
     85                         end
     86                         Dominated = m > M;
     87                         if Dominated || M == 2
     88                             break;
     89                         end
     90                     end
     91                 end
     92                 if ~Dominated
     93                     FrontNo(i) = MaxFNo;
     94                 end
     95             end
     96         end
     97     end
     98     FrontNo(rank) = FrontNo;
     99     FrontNo = FrontNo(:,Loc);
    100 end
    101 
    102 function [FrontNo,MaxFNo] = T_ENS(PopObj,nSort)
    103     [PopObj,~,Loc] = unique(PopObj,'rows');
    104     Table          = hist(Loc,1:max(Loc));
    105     [N,M]          = size(PopObj);
    106     [PopObj,rank]  = sortrows(PopObj);
    107     FrontNo        = inf(1,N);
    108     MaxFNo         = 0;
    109     Forest    = zeros(1,N);
    110     Children  = zeros(N,M-1);
    111     LeftChild = zeros(1,N) + M;
    112     Father    = zeros(1,N);
    113     Brother   = zeros(1,N) + M;
    114     [~,ORank] = sort(PopObj(:,2:M),2,'descend');
    115     ORank     = ORank + 1;
    116     while sum(Table(FrontNo<inf)) < min(nSort,length(Loc))
    117         MaxFNo = MaxFNo + 1;
    118         root   = find(FrontNo==inf,1);
    119         Forest(MaxFNo) = root;
    120         FrontNo(root)  = MaxFNo;
    121         for p = 1 : N
    122             if FrontNo(p) == inf
    123                 Pruning = zeros(1,N);
    124                 q = Forest(MaxFNo);
    125                 while true
    126                     m = 1;
    127                     while m < M && PopObj(p,ORank(q,m)) >= PopObj(q,ORank(q,m))
    128                         m = m + 1;
    129                     end
    130                     if m == M
    131                         break;
    132                     else
    133                         Pruning(q) = m;
    134                         if LeftChild(q) <= Pruning(q)
    135                             q = Children(q,LeftChild(q));
    136                         else
    137                             while Father(q) && Brother(q) > Pruning(Father(q))
    138                                 q = Father(q);
    139                             end
    140                             if Father(q)
    141                                 q = Children(Father(q),Brother(q));
    142                             else
    143                                 break;
    144                             end
    145                         end
    146                     end
    147                 end
    148                 if m < M
    149                     FrontNo(p) = MaxFNo;
    150                     q = Forest(MaxFNo);
    151                     while Children(q,Pruning(q))
    152                         q = Children(q,Pruning(q));
    153                     end
    154                     Children(q,Pruning(q)) = p;
    155                     Father(p) = q;
    156                     if LeftChild(q) > Pruning(q)
    157                         Brother(p)   = LeftChild(q);
    158                         LeftChild(q) = Pruning(q);
    159                     else
    160                         bro = Children(q,LeftChild(q));
    161                         while Brother(bro) < Pruning(q)
    162                             bro = Children(q,Brother(bro));
    163                         end
    164                         Brother(p)   = Brother(bro);
    165                         Brother(bro) = Pruning(q);
    166                     end
    167                 end
    168             end
    169         end
    170     end
    171     FrontNo(rank) = FrontNo;
    172     FrontNo = FrontNo(:,Loc);
    173 end

    IGD.m

     1 function Score = IGD(PopObj,PF)
     2 % <metric> <min>
     3 % Inverted generational distance
     4 
     5 %------------------------------- Reference --------------------------------
     6 % C. A. Coello Coello and N. C. Cortes, Solving multiobjective optimization
     7 % problems using an artificial immune system, Genetic Programming and
     8 % Evolvable Machines, 2005, 6(2): 163-190.
     9 %------------------------------- Copyright --------------------------------
    10 % Copyright (c) 2018-2019 BIMK Group. You are free to use the PlatEMO for
    11 % research purposes. All publications which use this platform or any code
    12 % in the platform should acknowledge the use of "PlatEMO" and reference "Ye
    13 % Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
    14 % for evolutionary multi-objective optimization [educational forum], IEEE
    15 % Computational Intelligence Magazine, 2017, 12(4): 73-87".
    16 %--------------------------------------------------------------------------
    17 
    18     Distance = min(pdist2(PF,PopObj),[],2);
    19     Score    = mean(Distance);
    20 end

    NSGAIII.m

     1 function NSGAIII(Global)
     2 
     3     %% Generate the reference points and random population
     4     [Z,Global.N] = UniformPoint(Global.N,Global.M);
     5     Population   = Global.Initialization();
     6     Zmin         = min(Population(all(Population.cons<=0,2)).objs,[],1);
     7 
     8     %% Optimization
     9     while Global.NotTermination(Population)
    10         MatingPool = TournamentSelection(2,Global.N,sum(max(0,Population.cons),2));
    11         Offspring  = GA(Population(MatingPool));
    12         Zmin       = min([Zmin;Offspring(all(Offspring.cons<=0,2)).objs],[],1);
    13         Population = EnvironmentalSelection([Population,Offspring],Global.N,Z,Zmin);
    14     end
    15 end

     

     

  • 相关阅读:
    Beetl模板 [记录]
    wx 小程序开发 [记录]
    高德定位获取省市区[记录]
    vue 学习记录 [记录]
    正则表达+验证 [记录]
    倒计时60s短信 [记录]
    @media [记录]
    JSON + Ajax [记录]
    Webstorm [记录]
    JQ 组合代码 [记录]
  • 原文地址:https://www.cnblogs.com/Vae1990Silence/p/12769336.html
Copyright © 2011-2022 走看看