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

 

 

原文地址:https://www.cnblogs.com/Vae1990Silence/p/12769336.html