• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

NSGA3理解(NSGA3算法及其MATLAB版本实现)

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

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  
                       
                    
                    

鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap