NSGA-Ⅲ源代码如下，供大家学习和应用。该算法在梯级水电-火电的应用订阅专栏即可查看：
1、主函数
%
% Copyright (c) 2016, Mostapha Kalami Heris & Yarpiz (www.yarpiz.com)
% All rights reserved. Please read the "LICENSE" file for license terms.
%
% Project Code: YPEA126
% Project Title: Non-dominated Sorting Genetic Algorithm III (NSGA-III)
% Publisher: Yarpiz (www.yarpiz.com)
%
% Implemented by: Mostapha Kalami Heris, PhD (member of Yarpiz Team)
%
% Cite as:
% Mostapha Kalami Heris, NSGA-III: Non-dominated Sorting Genetic Algorithm, the Third Version — MATLAB Implementation (URL: https://yarpiz.com/456/ypea126-nsga3), Yarpiz, 2016.
%
% Contact Info: sm.kalami@gmail.com, info@yarpiz.com
%
% Base Reference Paper:
% K. Deb and H. Jain, "An Evolutionary Many-Objective Optimization Algorithm
% Using Reference-Point-Based Nondominated Sorting Approach, Part I: Solving
% Problems With Box Constraints, "
% in IEEE Transactions on Evolutionary Computation,
% vol. 18, no. 4, pp. 577-601, Aug. 2014.
%
% Reference Paper URL: http://doi.org/10.1109/TEVC.2013.2281535
%
nsga3;
2、NSGA-Ⅲ函数
clc;
clear;
close all;
%% Problem Definition
CostFunction = @(x) MOP2(x); % Cost Function
nVar = 5; % Number of Decision Variables
VarSize = [1 nVar]; % Size of Decision Variables Matrix
VarMin = -1; % Lower Bound of Variables
VarMax = 1; % Upper Bound of Variables
% Number of Objective Functions
nObj = numel(CostFunction(unifrnd(VarMin, VarMax, VarSize)));
%% NSGA-II Parameters
% Generating Reference Points
nDivision = 10;
Zr = GenerateReferencePoints(nObj, nDivision);
MaxIt = 50; % Maximum Number of Iterations
nPop = 80; % Population Size
pCrossover = 0.5; % Crossover Percentage
nCrossover = 2*round(pCrossover*nPop/2); % Number of Parnets (Offsprings)
pMutation = 0.5; % Mutation Percentage
nMutation = round(pMutation*nPop); % Number of Mutants
mu = 0.02; % Mutation Rate
sigma = 0.1*(VarMax-VarMin); % Mutation Step Size
%% Colect Parameters
params.nPop = nPop;
params.Zr = Zr;
params.nZr = size(Zr, 2);
params.zmin = [];
params.zmax = [];
params.smin = [];
%% Initialization
disp('Staring NSGA-III ...');
empty_individual.Position = [];
empty_individual.Cost = [];
empty_individual.Rank = [];
empty_individual.DominationSet = [];
empty_individual.DominatedCount = [];
empty_individual.NormalizedCost = [];
empty_individual.AssociatedRef = [];
empty_individual.DistanceToAssociatedRef = [];
pop = repmat(empty_individual, nPop, 1);
for i = 1:nPop
pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
pop(i).Cost = CostFunction(pop(i).Position);
end
% Sort Population and Perform Selection
[pop, F, params] = SortAndSelectPopulation(pop, params);
%% NSGA-II Main Loop
for it = 1:MaxIt
% Crossover
popc = repmat(empty_individual, nCrossover/2, 2);
for k = 1:nCrossover/2
i1 = randi([1 nPop]);
p1 = pop(i1);
i2 = randi([1 nPop]);
p2 = pop(i2);
[popc(k, 1).Position, popc(k, 2).Position] = Crossover(p1.Position, p2.Position);
popc(k, 1).Cost = CostFunction(popc(k, 1).Position);
popc(k, 2).Cost = CostFunction(popc(k, 2).Position);
end
popc = popc(:);
% Mutation
popm = repmat(empty_individual, nMutation, 1);
for k = 1:nMutation
i = randi([1 nPop]);
p = pop(i);
popm(k).Position = Mutate(p.Position, mu, sigma);
popm(k).Cost = CostFunction(popm(k).Position);
end
% Merge
pop = [pop
popc
popm]; %#ok
% Sort Population and Perform Selection
[pop, F, params] = SortAndSelectPopulation(pop, params);
% Store F1
F1 = pop(F{1});
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Number of F1 Members = ' num2str(numel(F1))]);
% Plot F1 Costs
figure(1);
PlotCosts(F1);
pause(0.01);
end
%% Results
disp(['Final Iteration: Number of F1 Members = ' num2str(numel(F1))]);
disp('Optimization Terminated.');
3、其他函数
3.1 AssociateToReferencePoint
function [pop, d, rho] = AssociateToReferencePoint(pop, params)
Zr = params.Zr;
nZr = params.nZr;
rho = zeros(1, nZr);
d = zeros(numel(pop), nZr);
for i = 1:numel(pop)
for j = 1:nZr
w = Zr(:, j)/norm(Zr(:, j));
z = pop(i).NormalizedCost;
d(i, j) = norm(z - w'*z*w);
end
[dmin, jmin] = min(d(i, :));
pop(i).AssociatedRef = jmin;
pop(i).DistanceToAssociatedRef = dmin;
rho(jmin) = rho(jmin) + 1;
end
end
3.2 Crossover
function [y1 y2] = Crossover(x1, x2)
alpha = rand(size(x1));
y1 = alpha.*x1+(1-alpha).*x2;
y2 = alpha.*x2+(1-alpha).*x1;
end
3.3 Dominate
function b = Dominates(x, y)
if isstruct(x)
x = x.Cost;
end
if isstruct(y)
y = y.Cost;
end
b = all(x <= y) && any(x<y);
end
3.4 GenerateReferencePoints
function Zr = GenerateReferencePoints(M, p)
Zr = GetFixedRowSumIntegerMatrix(M, p)' / p;
end
function A = GetFixedRowSumIntegerMatrix(M, RowSum)
if M < 1
error('M cannot be less than 1.');
end
if floor(M) ~= M
error('M must be an integer.');
end
if M == 1
A = RowSum;
return;
end
A = [];
for i = 0:RowSum
B = GetFixedRowSumIntegerMatrix(M - 1, RowSum - i);
A = [A; i*ones(size(B, 1), 1) B];
end
end
3.5 MOP2
function z = MOP2(x)
n = numel(x);
z1 = 1-exp(-sum((x-1/sqrt(n)).^2));
z2 = 1-exp(-sum((x+1/sqrt(n)).^2));
z = [z1 z2]';
end
3.6 Mutate
function y = Mutate(x, mu, sigma)
nVar = numel(x);
nMu = ceil(mu*nVar);
j = randsample(nVar, nMu);
y = x;
y(j) = x(j)+sigma*randn(size(j));
end
3.7 NonDominatedSorting
function [pop, F] = NonDominatedSorting(pop)
nPop = numel(pop);
for i = 1:nPop
pop(i).DominationSet = [];
pop(i).DominatedCount = 0;
end
F{1} = [];
for i = 1:nPop
for j = i+1:nPop
p = pop(i);
q = pop(j);
if Dominates(p, q)
p.DominationSet = [p.DominationSet j];
q.DominatedCount = q.DominatedCount+1;
end
if Dominates(q.Cost, p.Cost)
q.DominationSet = [q.DominationSet i];
p.DominatedCount = p.DominatedCount+1;
end
pop(i) = p;
pop(j) = q;
end
if pop(i).DominatedCount == 0
F{1} = [F{1} i];
pop(i).Rank = 1;
end
end
k = 1;
while true
Q = [];
for i = F{k}
p = pop(i);
for j = p.DominationSet
q = pop(j);
q.DominatedCount = q.DominatedCount-1;
if q.DominatedCount == 0
Q = [Q j];
q.Rank = k+1;
end
pop(j) = q;
end
end
if isempty(Q)
break;
end
F{k+1} = Q;
k = k+1;
end
end
3.8 NormalizePopulation
function [pop, params] = NormalizePopulation(pop, params)
params.zmin = UpdateIdealPoint(pop, params.zmin);
fp = [pop.Cost] - repmat(params.zmin, 1, numel(pop));
params = PerformScalarizing(fp, params);
a = FindHyperplaneIntercepts(params.zmax);
for i = 1:numel(pop)
pop(i).NormalizedCost = fp(:, i)./a;
end
end
function a = FindHyperplaneIntercepts(zmax)
w = ones(1, size(zmax, 2))/zmax;
a = (1./w)';
end
3.9 PerformScalarizing
function params = PerformScalarizing(z, params)
nObj = size(z, 1);
nPop = size(z, 2);
if ~isempty(params.smin)
zmax = params.zmax;
smin = params.smin;
else
zmax = zeros(nObj, nObj);
smin = inf(1, nObj);
end
for j = 1:nObj
w = GetScalarizingVector(nObj, j);
s = zeros(1, nPop);
for i = 1:nPop
s(i) = max(z(:, i)./w);
end
[sminj, ind] = min(s);
if sminj < smin(j)
zmax(:, j) = z(:, ind);
smin(j) = sminj;
end
end
params.zmax = zmax;
params.smin = smin;
end
function w = GetScalarizingVector(nObj, j)
epsilon = 1e-10;
w = epsilon*ones(nObj, 1);
w(j) = 1;
end
3.10 SortAndSelectPopulation
function [pop, F, params] = SortAndSelectPopulation(pop, params)
[pop, params] = NormalizePopulation(pop, params);
[pop, F] = NonDominatedSorting(pop);
nPop = params.nPop;
if numel(pop) == nPop
return;
end
[pop, d, rho] = AssociateToReferencePoint(pop, params);
newpop = [];
for l = 1:numel(F)
if numel(newpop) + numel(F{l}) > nPop
LastFront = F{l};
break;
end
newpop = [newpop; pop(F{l})]; %#ok
end
while true
[~, j] = min(rho);
AssocitedFromLastFront = [];
for i = LastFront
if pop(i).AssociatedRef == j
AssocitedFromLastFront = [AssocitedFromLastFront i]; %#ok
end
end
if isempty(AssocitedFromLastFront)
rho(j) = inf;
continue;
end
if rho(j) == 0
ddj = d(AssocitedFromLastFront, j);
[~, new_member_ind] = min(ddj);
else
new_member_ind = randi(numel(AssocitedFromLastFront));
end
MemberToAdd = AssocitedFromLastFront(new_member_ind);
LastFront(LastFront == MemberToAdd) = [];
newpop = [newpop; pop(MemberToAdd)]; %#ok
rho(j) = rho(j) + 1;
if numel(newpop) >= nPop
break;
end
end
[pop, F] = NonDominatedSorting(newpop);
end
3.11 UpdateIdealPoint
function zmin = UpdateIdealPoint(pop, prev_zmin)
if ~exist('prev_zmin', 'var') || isempty(prev_zmin)
prev_zmin = inf(size(pop(1).Cost));
end
zmin = prev_zmin;
for i = 1:numel(pop)
zmin = min(zmin, pop(i).Cost);
end
end
