小编在查阅关于NSGA-II相关的源代码时，在matlab官网时发现一个质量较高的源代码，点击左下角阅读原文可直接查看此源代码。如果各位想重新复习一下关于NSGA-II的理论知识，可以查看NSGA-II多目标优化算法讲解（附MATLAB代码）这篇推文。
这份源代码一共对14个测试函数进行测试，下面小编将这份源代码中的各个函数逐个列出来。
01 | Main_NSGA2
主函数Main_NSGA2，运行主函数的时候，命令行窗口会出现Test problem index :，这时需要输入1~14中的任意一个数字，意思就是选择14个测试函数中的任意一个函数。
%% https://ww2.mathworks.cn/matlabcentral/fileexchange/49806-matlab-code-for-constrained-nsga-ii-dr-s-baskar-s-tamilselvi-and-p-r-varshini
clear all
clc
global V M xl xu etac etam p pop_size pm
%% Description
% 1. This is the main program of NSGA II. It requires only one input, which is test problem
% index, 'p'. NSGA II code is tested and verified for 14 test problems.
% 2. This code defines population size in 'pop_size', number of design
% variables in 'V', number of runs in 'no_runs', maximum number of
% generations in 'gen_max', current generation in 'gen_count' and number of objectives
% in 'M'.
% 3. 'xl' and 'xu' are the lower and upper bounds of the design variables.
% 4. Final optimal Pareto soutions are in the variable 'pareto_rank1', with design
% variables in the coumns (1:V), objectives in the columns (V+1 to V+M),
% constraint violation in the column (V+M+1), Rank in (V+M+2), Distance in (V+M+3).
%% code starts
M=2;
p=input('Test problem index :');
pop_size=200; % Population size
no_runs=1; % Number of runs
gen_max=500; % MAx number of generations - stopping criteria
fname='test_case'; % Objective function and constraint eval（p==2 || p==5 || p==7), gen_max=1000; end
if p<=9 % Unconstrained test functions
tV=[2;30;3;1;30;4;30;10;10];
V=tV(p);
txl=[-5*ones(1,V);zeros(1,V);-5*ones(1,V);-1000*ones(1,V);zeros(1,V);-1/sqrt(V)*ones(1,V);zeros(1,V); 0 -5*ones(1,V-1);zeros(1,V)];
txu=[10*ones(1,V); ones(1,V);5*ones(1,V);1000*ones(1,V);ones(1,V);1/sqrt(V) *ones(1,V);ones(1,V);1 5*ones(1,V-1);ones(1,V)];
xl=(txl(p,1:V)); % lower bound vector
xu=(txu(p,1:V)); % upper bound vectorfor
etac = 20; % distribution index for crossover
etam = 20; % distribution index for mutation / mutation constant
else % Constrained test functions
p1=p-9;
tV=[2;2;2;6;2];
V=tV(p1);
txl=[0 0 0 0 0 0;-20 -20 0 0 0 0;0 0 0 0 0 0;0 0 1 0 1 0;0.1 0 0 0 0 0];
txu=[5 3 0 0 0 0;20 20 0 0 0 0;pi pi 0 0 0 0;10 10 5 6 5 10;1 5 0 0 0 0];
xl=(txl(p1,1:V)); % lower bound vector
xu=(txu(p1,1:V)); % upper bound vectorfor i=1:NN
etac = 20; % distribution index for crossover
etam = 100; % distribution index for mutation / mutation constant
end
pm=1/V; % Mutation Probability
Q=[];
for run = 1:no_runs
%% Initial population
xl_temp=repmat(xl, pop_size,1);
xu_temp=repmat(xu, pop_size,1);
x = xl_temp+((xu_temp-xl_temp).*rand(pop_size,V));
%% Evaluate objective function
for i =1:pop_size
[ff(i,:),err(i,:)] =feval（fname, x(i,:)); % Objective function evaulation
end
error_norm=normalisation(err); % Normalisation of the constraint violation
population_init=[x ff error_norm];
[population,front]=NDS_CD_cons(population_init); % Non domination Sorting on initial population
%% Generation Starts
for gen_count=1:gen_max
% selection (Parent Pt of 'N' pop size)
parent_selected=tour_selection(population); % 10 Tournament selection
%% Reproduction (Offspring Qt of 'N' pop size)
child_offspring = genetic_operator(parent_selected(:,1:V)); % SBX crossover and polynomial mutation
for ii = 1:pop_size
[fff(ii,:) err(ii,:)]=feval（fname, child_offspring(ii,:)); % objective function eval（err);
child_offspring=[child_offspring fff error_norm];
%% INtermediate population (Rt= Pt U Qt of 2N size)
population_inter=[population(:,1:V+M+1) ; child_offspring(:,1:V+M+1)];
[population_inter_sorted, front]=NDS_CD_cons(population_inter); % Non domination Sorting on offspring
%% Replacement - N
new_pop=replacement(population_inter_sorted, front);
population=new_pop;
end
new_pop=sortrows(new_pop,V+1);
paretoset(run).trial=new_pop(:,1:V+M+1);
Q = [Q; paretoset(run).trial]; % Combining Pareto solutions obtained in each run
end
%% Result and Pareto plot
if run==1
plot(new_pop(:,V+1),new_pop(:,V+2),'*')
else
[pareto_filter front]=NDS_CD_cons(Q); % Applying non domination sorting on the combined Pareto solution set
rank1_index=find(pareto_filter(:,V+M+2)==1); % Filtering the best solutions of rank 1 Pareto
pareto_rank1=pareto_filter(rank1_index,1:V+M)
plot(pareto_rank1(:,V+1),pareto_rank1(:,V+2),'*') % Final Pareto plot
end
xlabel('objective function 1')
ylabel('objective function 2')
if p==1
title(' 1 - Test case 1')
elseif p==2
title(' 2 - ZDT1')
elseif p==3
title(' 3 - KUR')
elseif p==4
title(' 4 - SCH')
elseif p==5
title(' 5 - ZDT2')
elseif p==6
title(' 6 - Test case 3')
elseif p==7
title(' 7 - ZDT3')
elseif p==8
title(' 8 - ZDT4')
elseif p==9
title(' 9 - ZDT6')
elseif p==10
title(' 10 - BNH')
elseif p==11
title(' 11 - SRN')
elseif p==12
title(' 12 - TNK')
elseif p==13
title(' 13 - OSY')
elseif p==14
title(' 14 - CONSTR')
end
02 | test_case
前9个测试函数是属于无约束优化函数，后5个测试函数属于有约束优化函数。返回值fit是目标函数值，返回值err是违反约束（超过约束）的值。比如说，一个函数的约束为x1+x2<=5，此时err=x1+x2-5，即当x1=4，x2=5时，err=4.
%% Description
% 1. This function returns the objective functions f1, and f2 in the vector 'fit' and
% constraints in the vector 'c' for the chromosome 'x'.
% 2. 'V' is the number of optimization variables.
% 3. All the constraints 'c' are converted to the form h(x)<=0.
% 4. Nine unconstrained test probems are used (p=1 to p=9)
% 5. Five constrained test problems are used (p=10 to p=14)
% 6. Refer above references for the details of each test problem: number of objectives, number of design variabes, their lower and upper limits,
% number of constraints, type of constraints etc,.
%% reference
% 1. BINH, Thanh. "A multiobjective evolutionary algorithm. The study cases".
% Technical report. Barleben, Germany. 1999.
% 2. DEB, Kalyanmoy. "Multi-Objective optimization using evolutionary
% algorithms". John Wiley & Sons, LTD. Kanpur, India. 2004.
function [fit err]=test_case(x)
global p V
%% Unconstrained Test functions (for p=1 to p=9)
if p==1 % Test case problem 1
f1=(4*x(1)^2)+(4*x(2)^2);
f2=((x(1)-5)^2)+((x(2)-5)^2);
end
if p==2 % ZDT1 from Deb paper NSGA2
cons=[0];
f1 = x(1);
g=1+(9*sum(x(2:V),2)/(V-1));
f2 = g*(1-sqrt(x(1)/g));
end
if p==3 % kUR from Deb
f1=(-10*exp(-0.2*(sqrt(x(1)^2+x(2)^2))))+(-10*exp(-0.2*(sqrt(x(2)^2+x(3)^2))));
f2=((abs(x(1))^0.8) + (5*sin(x(1))^3))+((abs(x(2))^0.8) + (5*sin(x(2))^3))+((abs(x(3))^0.8) + (5*sin(x(3))^3));
end
if p==4 % SCH frm Deb paper
f1=x.*x;
f2=(x-2).^2;
end
if p==5 % ZDT2
f1 = (x(1));
g=1+(9*sum(x(2:V),2)/(V-1));
f2 =((1-(x(1)/g)^2));
end
if p==6 % Test case problem 2
f1=1-exp(-sum((x-1/sqrt(V)).^2,2));
f2=1-exp(-sum((x+1/sqrt(V)).^2,2));
end
if p==7 % ZDT3
f1 = x(1);
g=1+(9*sum(x(2:V),2)/(V-1));
f2 = (1-(sqrt(x(1)/g)) - ((x(1)/g)*sin(10*pi*x(1))));
end
if p==8 % ZDT4
f1 = x(1); temp=0;
for ii = 2: V
temp=temp+((x(ii))^2)-(10*cos(4*pi*x(ii)));
end
g= 1 + (10*(V-1)) + temp;
f2 = (1-sqrt(x(1)/g));
end
if p==9 % ZDT6
f1 = 1-(exp(-4*x(1)))*(sin(6*pi*x(1)))^6;
g=1+(9*(sum(x(2:V),2)/(V-1))^0.25);
f2 = (1-(f1/g)^2);
end
err= zeros(1,1);
%% Constrained Test functions (for p=10 to p=14)
if p==10 %BNH
f1=4*(x(1)^2)+4*(x(2)^2);
f2=(x(1)-5)^2+(x(2)-5)^2;
c(1,1)=(x(1)-5)^2 + x(2)^2 -25;
c(1,2)=-(x(1)-8)^2-(x(2)+3)^2+7.7;
err=(c>0).*c;
end
if p==11 %SRN
f1=(x(1)-2)^2+(x(2)-1)^2+2;
f2=9*x(1)-(x(2)-1)^2;
c(1,1)=x(1)^2+x(2)^2-225;
c(1,2)=x(1)-(3*x(2))+10;
err=(c>0).*c;
end
if p==12 %TNK
f1=x(1);
f2=x(2);
c(1,1)=-x(1)^2-x(2)^2+1+(0.1*cos(16*atan((x(1)/x(2)))));
c(1,2)=(x(1)-0.5)^2+(x(2)-0.5)^2-0.5;
err=(c>0).*c;
end
if p==13 % OSY
f1=-((25*(x(1)-2)^2)+((x(2)-2)^2)+((x(3)-1)^2)+((x(4)-4)^2)+((x(5)-1)^2));
f2=(x(1)^2)+(x(2)^2)+(x(3)^2)+(x(4)^2)+(x(5)^2)+(x(6)^2);
c(1,1)=-x(1)-x(2)+2;
c(1,2)=-6+x(1)+x(2);
c(1,3)=-2+x(2)-x(1);
c(1,4)=-2+x(1)-3*x(2);
c(1,5)=-4+((x(3)-3)^2)+x(4);
c(1,6)=-((x(5)-3)^2)-x(6)+4;
err=(c>0).*c;
end
if p==14 % CONSTR
f1=x(1);
f2=(1+x(2))/(x(1));
c(1,1)=-x(2)-(9*x(1))+6;
c(1,2)=+x(2)-9*x(1)+1;
err=(c>0).*c;
end
fit=[f1 f2];
03 | normalisation
normalisation函数将违反约束的数值标准化。比如说输入[1,4;2,3]，那么这两列中每一列最大值的矩阵为[2,4]；然后，将第第一列中每个数除以2，第二列中每个数除以4，则标准化成[0.5,1;1,0.75]；最后，将每一行的数值相加，则变成[1.5;1.75]。
function err_norm = normalisation(error_pop)
%% Description
% 1. This function normalises the constraint violation of various individuals, since the range of
% constraint violation of every chromosome is not uniform.
% 2. Input is in the matrix error_pop with size [pop_size, number of constraints].
% 3. Output is a normalised vector, err_norm of size [pop_size,1]
%% Error Nomalisation
[N,nc]=size(error_pop);
con_max=0.001+max(error_pop);
con_maxx=repmat(con_max,N,1);
cc=error_pop./con_maxx;
err_norm=sum(cc,2); % finally sum up all violations
04 | NDS_CD_cons
这是关于非支配排序的函数，具体的执行步骤同NSGA-II多目标优化算法讲解（附MATLAB代码）这篇推文的执行步骤。
%% Description
% 1. This function is to perform Deb's fast elitist non-domination sorting and crowding distance assignment.
% 2. Input is in the variable 'population' with size: [size(popuation), V+M+1]
% 3. This function returns 'chromosome_NDS_CD' with size [size(population),V+M+3]
% 4. A flag 'problem_type' is used to identify whether the population is fully feasible (problem_type=0) or fully infeasible (problem_type=1)
% or partly feasible (problem_type=0.5).
%% Reference:
%Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and T. Meyarivan, " A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II",
%IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 6, No. 2, APRIL 2002.
%% function begins
function [chromosome_NDS_CD front] = NDS_CD_cons(population)
global V M problem_type
%% Initialising structures and variables
chromosome_NDS_CD1=[];
infpop=[];
front.fr=[];
struct.sp=[];
rank=1;
%% Segregating feasible and infeasible solutions
if all(population(:,V+M+1)==0)
problem_type=0;
chromosome=population(:,1:V+M); % All Feasible chromosomes;
pop_size1=size(chromosome,1);
elseif all(population(:,V+M+1)~=0)
problem_type=1;
pop_size1=0;
infchromosome=population; % All InFeasible chromosomes;
else
problem_type=0.5;
feas_index=find(population(:,V+M+1)==0);
chromosome=population(feas_index,1:V+M); % Feasible chromosomes;
pop_size1=size(chromosome,1);
infeas_index=find(population(:,V+M+1)~=0);
infchromosome=population(infeas_index,1:V+M+1); % infeasible chromosomes;
end
%% Handling feasible solutions
if problem_type==0 | problem_type==0.5
pop_size1 = size(chromosome,1);
f1 = chromosome(:,V+1); % objective function values
f2 = chromosome(:,V+2);
%Non- Domination Sorting
% First front
for p=1:pop_size1
struct(p).sp=find(((f1(p)-f1)<0 &(f2(p)-f2)<0) | ((f2(p)-f2)==0 &(f1(p)-f1)<0) | ((f1(p)-f1)==0 &(f2(p)-f2)<0));
n(p)=length(find(((f1(p)-f1)>0 &(f2(p)-f2)>0) | ((f2(p)-f2)==0 &(f1(p)-f1)>0) | ((f1(p)-f1)==0 &(f2(p)-f2)>0)));
end
front(1).fr=find(n==0);
% Creating subsequent fronts
while (~isempty(front(rank).fr))
front_indiv=front(rank).fr;
n(front_indiv)=inf;
chromosome(front_indiv,V+M+1)=rank;
rank=rank+1;
front(rank).fr=[];
for i = 1:length(front_indiv)
temp=struct(front_indiv(i)).sp;
n(temp)=n(temp)-1;
end
q=find(n==0);
front(rank).fr=[front(rank).fr q];
end
chromosome_sorted=sortrows(chromosome,V+M+1); % Ranked population
%Crowding distance Assignment
rowsindex=1;
for i = 1:length(front)-1
l_f=length(front(i).fr);
if l_f > 2
sorted_indf1=[];
sorted_indf2=[];
sortedf1=[];
sortedf2=[];
% sorting based on f1 and f2;
[sortedf1 sorted_indf1]=sortrows(chromosome_sorted(rowsindex:(rowsindex+l_f-1),V+1));
[sortedf2 sorted_indf2]=sortrows(chromosome_sorted(rowsindex:(rowsindex+l_f-1),V+2));
f1min=chromosome_sorted(sorted_indf1(1)+rowsindex-1,V+1);
f1max=chromosome_sorted(sorted_indf1(end)+rowsindex-1,V+1);
chromosome_sorted(sorted_indf1(1)+rowsindex-1,V+M+2)=inf;
chromosome_sorted(sorted_indf1(end)+rowsindex-1,V+M+2)=inf;
f2min=chromosome_sorted(sorted_indf2(1)+rowsindex-1,V+2);
f2max=chromosome_sorted(sorted_indf2(end)+rowsindex-1,V+2);
chromosome_sorted(sorted_indf2(1)+rowsindex-1,V+M+3)=inf;
chromosome_sorted(sorted_indf2(end)+rowsindex-1,V+M+3)=inf;
for j = 2:length(front(i).fr)-1
if (f1max - f1min == 0) | (f2max - f2min == 0)
chromosome_sorted(sorted_indf1(j)+rowsindex-1,V+M+2)=inf;
chromosome_sorted(sorted_indf2(j)+rowsindex-1,V+M+3)=inf;
else
chromosome_sorted(sorted_indf1(j)+rowsindex-1,V+M+2)=(chromosome_sorted(sorted_indf1(j+1)+rowsindex-1,V+1)-chromosome_sorted(sorted_indf1(j-1)+rowsindex-1,V+1))/(f1max-f1min);
chromosome_sorted(sorted_indf2(j)+rowsindex-1,V+M+3)=(chromosome_sorted(sorted_indf2(j+1)+rowsindex-1,V+2)-chromosome_sorted(sorted_indf2(j-1)+rowsindex-1,V+2))/(f2max-f2min);
end
end
else
chromosome_sorted(rowsindex:(rowsindex+l_f-1),V+M+2:V+M+3)=inf;
end
rowsindex = rowsindex + l_f;
end
chromosome_sorted(:,V+M+4) = sum(chromosome_sorted(:,V+M+2:V+M+3),2);
chromosome_NDS_CD1 = [chromosome_sorted(:,1:V+M) zeros(pop_size1,1) chromosome_sorted(:,V+M+1) chromosome_sorted(:,V+M+4)]; % Final Output Variable
end
%% Handling infeasible solutions
if problem_type==1 | problem_type==0.5
infpop=sortrows(infchromosome,V+M+1);
infpop=[infpop(:,1:V+M+1) (rank:rank-1+size(infpop,1))' inf*(ones(size(infpop,1),1))];
for kk = (size(front,2)):(size(front,2))+(length(infchromosome))-1;
front(kk).fr= pop_size1+1;
end
end
chromosome_NDS_CD = [chromosome_NDS_CD1;infpop];
05 | tour_selection
function [parent_selected] = tour_selection(pool)
%% Description
% 1. Parents are selected from the population pool for reproduction by using binary tournament selection
% based on the rank and crowding distance.
% 2. An individual is selected if the rank is lesser than the other or if
% crowding distance is greater than the other.
% 3. Input and output are of same size [pop_size, V+M+3].
%% Binary Tournament Selection
[pop_size, distance]=size(pool);
rank=distance-1;
candidate=[randperm(pop_size);randperm(pop_size)]';
for i = 1: pop_size
parent=candidate(i,:); % Two parents indexes are randomly selected
if pool(parent(1),rank)~=pool(parent(2),rank) % For parents with different ranks
if pool(parent(1),rank)<pool(parent(2),rank) % Checking the rank of two individuals
mincandidate=pool(parent(1),:);
elseif pool(parent(1),rank)>pool(parent(2),rank)
mincandidate=pool(parent(2),:);
end
parent_selected(i,:)=mincandidate; % Minimum rank individual is selected finally
else % for parents with same ranks
if pool(parent(1),distance)>pool(parent(2),distance) % Checking the distance of two parents
maxcandidate=pool(parent(1),:);
elseif pool(parent(1),distance)< pool(parent(2),distance)
maxcandidate=pool(parent(2),:);
else
temp=randperm(2);
maxcandidate=pool(parent(temp(1)),:);
end
parent_selected(i,:)=maxcandidate; % Maximum distance individual is selected finally
end
end
06 | genetic_operator
function child_offspring = genetic_operator(parent_selected)
global V xl xu etac
%% Description
% 1. Crossover followed by mutation
% 2. Input is in 'parent_selected' matrix of size [pop_size,V].
% 3. Output is also of same size in 'child_offspring'.
%% Reference
% Deb & samir agrawal,"A Niched-Penalty Approach for Constraint Handling in Genetic Algorithms".
%% SBX cross over operation incorporating boundary constraint
[N] = size(parent_selected,1);
xl1=xl';
xu1=xu';
rc=randperm(N);
for i=1:(N/2)
parent1=parent_selected((rc(2*i-1)),:);
parent2=parent_selected((rc(2*i)),:);
if (isequal(parent1,parent2))==1 & rand(1)>0.5
child1=parent1;
child2=parent2;
else
for j = 1: V
if parent1(j)<parent2(j)
beta(j)= 1 + (2/(parent2(j)-parent1(j)))*(min((parent1(j)-xl1(j)),(xu1(j)-parent2(j))));
else
beta(j)= 1 + (2/(parent1(j)-parent2(j)))*(min((parent2(j)-xl1(j)),(xu1(j)-parent1(j))));
end
end
u=rand(1,V);
alpha=2-beta.^-(etac+1);
betaq=(u<=(1./alpha)).*(u.*alpha).^(1/(etac+1))+(u>(1./alpha)).*(1./(2 - u.*alpha)).^(1/(etac+1));
child1=0.5*(((1).*parent1) + (1).*parent2);
child2=0.5*(((1).*parent1) + (1).*parent2);
end
child_offspring((rc(2*i-1)),:)=poly_mutation(child1); % polynomial mutation
child_offspring((rc(2*i)),:)=poly_mutation(child2); % polynomial mutation
end
07 | poly_mutation
function mutated_child = poly_mutation(y)
global V xl xu etam pm
%% Description
% 1. Input is the crossovered child of size (1,V) in the vector 'y' from 'genetic_operator.m'.
% 2. Output is in the vector 'mutated_child' of size (1,V).
%% Polynomial mutation including boundary constraint
del=min((y-xl),(xu-y))./(xu-xl);
t=rand(1,V);
loc_mut=t<pm;
u=rand(1,V);
delq=(u<=0.5).*((((2*u)+((1-2*u).*((1-del).^(etam+1)))).^(1/(etam+1)))-1)+(u>0.5).*(1-((2*(1-u))+(2*(u-0.5).*((1-del).^(etam+1)))).^(1/(etam+1)));
c=y+delq.*loc_mut.*(xu-xl);
mutated_child=c;
08 | replacement
function new_pop=replacement(population_inter_sorted, front)
global pop_size
%% Description
% The next generation population is formed by appending each front subsequently until the
% population size exceeds the current population size. If When adding all the individuals
% of any front, the population exceeds the population size, then the required number of
% remaining individuals alone are selected from that particular front based
% on crowding distance.
%% code starts
index=0;
ii=1;
while index < pop_size
l_f=length(front(ii).fr);
if index+l_f < pop_size
new_pop(index+1:index+l_f,:)= population_inter_sorted(index+1:index+l_f,:);
index=index+l_f;
else
temp1=population_inter_sorted(index+1:index+l_f,:);
temp2=sortrows(temp1,size(temp1,2));
new_pop(index+1:pop_size,:)= temp2(l_f-(pop_size-index)+1:l_f,:);
index=index+l_f;
end
ii=ii+1;
end
09 | reference
1. Tamilselvi Selvaraj (2020). MATLAB code for Constrained NSGA II - Dr.S.Baskar, S. Tamilselvi and P.R.Varshini (https://www.mathworks.com/matlabcentral/fileexchange/49806-matlab-code-for-constrained-nsga-ii-dr-s-baskar-s-tamilselvi-and-p-r-varshini), MATLAB Central File Exchange. Retrieved January 20, 2020.
微博/知乎：随心390
长按识别二维码关注我们
