NSGA-ⅡMATLAB代码（转载）

function f = crowding_distance(x,problem)

% This function calculates the crowding distance

[N,M] = size(x);

switch problem

    case 1

       M = 2;

       V = 6;

    case 2

       M = 3;

       V = 12;

end

% Crowding distance for each front

for i = 1 : length(F(front).f)

    y(i,:) = x(F(front).f(i),:);

end

for i = 1 : M

    [sorted(i).individual,sorted(i).index] = sort(y(:,V + i));

    distance(sorted(i).index(1)).individual = Inf;

    distance(sorted(i).index(length(sorted(i).index))).individual = Inf;

end

[num,len] = size(y);

% Initialize all the distance of individuals as zero.

for i = 1 : M

    for j = 2 : num - 1

       distance(j).individual = 0;

    end

    objective(i).range = ...

               sorted(i).individual(length(sorted(i).individual)) - ...

               sorted(i).individual(1);

       % Maximum and minimum objectives value for the ith objective

end

% Caluclate the crowding distance for front one.

for i = 1 : M

    for j = 2 : num - 1

       distance(j).individual = distance(j).individual + ...

           (sorted(i).individual(j + 1) - sorted(i).individual(j - 1))/...

           objective(i).range;

       y(sorted(i).index(j),M + V + 2) = distance(j).individual;

    end

end

</PRE

        Published withMATLAB® 7.0

function f =genetic_operator(parent_chromosome,pro,mu,mum);

% This function is utilized to produce offsprings from parent chromosomes.

% The genetic operators corssover and mutation which are carried out with

% slight modifications from the original design. For more information read

% the document enclosed.

[N,M] = size(parent_chromosome);

switch pro

    case 1

       M = 2;

       V = 6;

    case 2

       M = 3;

       V = 12;

end

p = 1;

was_crossover = 0;

was_mutation = 0;

l_limit = 0;

u_limit = 1;

for i = 1 : N

    if rand(1) < 0.9

       child_1 = [];

       child_2 = [];

       parent_1 = round(N*rand(1));

       if parent_1 < 1

           parent_1 = 1;

       end

       parent_2 = round(N*rand(1));

       if parent_2 < 1

           parent_2 = 1;

       end

       while isequal(parent_chromosome(parent_1,:),parent_chromosome(parent_2,:))

           parent_2 = round(N*rand(1));

           if parent_2 < 1

               parent_2 = 1;

           end

       end

       parent_1 = parent_chromosome(parent_1,:);

       parent_2 = parent_chromosome(parent_2,:);

       for j = 1 : V

           % SBX (Simulated Binary Crossover)

           % Generate a random number

           u(j) = rand(1);

           if u(j) <= 0.5

               bq(j) = (2*u(j))^(1/(mu+1));

           else

               bq(j) = (1/(2*(1 - u(j))))^(1/(mu+1));

           end

           child_1(j) = ...

               0.5*(((1 + bq(j))*parent_1(j)) + (1 - bq(j))*parent_2(j));

           child_2(j) = ...

               0.5*(((1 - bq(j))*parent_1(j)) + (1 + bq(j))*parent_2(j));

           if child_1(j) > u_limit

               child_1(j) = u_limit;

           elseif child_1(j) < l_limit

               child_1(j) = l_limit;

           end

           if child_2(j) > u_limit

               child_2(j) = u_limit;

           elseif child_2(j) < l_limit

               child_2(j) = l_limit;

           end

       end

       child_1(:,V + 1: M + V) = evaluate_objective(child_1,pro);

       chi
4000
ld_2(:,V + 1: M + V) = evaluate_objective(child_2,pro);

       was_crossover = 1;

       was_mutation = 0;

    else

       parent_3 = round(N*rand(1));

       if parent_3 < 1

           parent_3 = 1;

       end

       % Make sure that the mutation does not result in variables out of

       % the search space. For both the MOP's the range for decision space

       % is [0,1]. In case different variables have different decision

       % space each variable can be assigned a range.

       child_3 = parent_chromosome(parent_3,:);

       for j = 1 : V

          r(j) = rand(1);

          if r(j) < 0.5

              delta(j) = (2*r(j))^(1/(mum+1)) - 1;

          else

              delta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1));

          end

          child_3(j) = child_3(j) + delta(j);

          if child_3(j) > u_limit

              child_3(j) = u_limit;

          elseif child_3(j) < l_limit

              child_3(j) = l_limit;

          end

       end

       child_3(:,V + 1: M + V) = evaluate_objective(child_3,pro);

        was_mutation = 1;

       was_crossover = 0;

    end

    if was_crossover

       child(p,:) = child_1;

       child(p+1,:) = child_2;

       was_cossover = 0;

       p = p + 2;

    elseif was_mutation

       child(p,:) = child_3(1,1 : M + V);

        was_mutation = 0;

       p = p + 1;

    end

end

f = child;

Published with MATLAB® 7.0

function f = initialize_variables(N,problem)

% function f = initialize_variables(N,problem)

% N - Population size

% problem - takes integer values 1 and 2 where,

%           '1' for MOP1

%          '2' for MOP2

%

% This function initializes the population with N individuals and each

% individual having M decision variables based on the selected problem.

% M = 6 for problem MOP1 and M = 12 for problem MOP2. The objective space

% for MOP1 is 2 dimensional while for MOP2 is 3 dimensional.

% Both the MOP's has 0 to 1 as its range for all the decision variables.

min = 0;

max = 1;

switch problem

    case 1

       M = 6;

       K = 8;

    case 2

       M = 12;

       K = 15;

end

for i = 1 : N

    % Initialize the decision variables

    for j = 1 : M

       f(i,j) = rand(1); % i.e f(i,j) = min + (max - min)*rand(1);

    end

    % Evaluate the objective function

    f(i,M + 1: K) = evaluate_objective(f(i,:),problem);

end

</PRE

        Published withMATLAB® 7.0

Non-Donimation Sort

This function sort the current popultion based onnon-domination. All the individuals in the first front are given a rank of 1,the second front individuals are assigned rank 2 and so on. After assigning therank the crowding in each front is calculated.

[N,M] = size(x);

switch problem

    case 1

       M = 2;

       V = 6;

    case 2

       M = 3;

       V = 12;

end

front = 1;

% There is nothing to this assignment, used only to manipulate easily in

% MATLAB.

F(front).f = [];

individual = [];

for i = 1 : N

    % Number of individuals that dominate this individual

    individual(i).n = 0;

    % Individuals which this individual dominate

    individual(i).p = [];

    for j = 1 : N

       dom_less = 0;

       dom_equal = 0;

       dom_more = 0;

       for k = 1 : M

           if (x(i,V + k) < x(j,V + k))

               dom_less = dom_less + 1;

           elseif (x(i,V + k) == x(j,V + k))

               dom_equal = dom_equal + 1;

           else

               dom_more = dom_more + 1;

           end

       end

       if dom_less == 0 & dom_equal ~= M

           individual(i).n = individual(i).n + 1;

       elseif dom_more == 0 & dom_equal ~= M

           individual(i).p = [individual(i).p j];

       end

    end

    if individual(i).n == 0

       x(i,M + V + 1) = 1;

       F(front).f = [F(front).f i];

    end

end

% Find the subsequent fronts

while ~isempty(F(front).f)

   Q = [];

   for i = 1 : length(F(front).f)

       if ~isempty(individual(F(front).f(i)).p)

           for j = 1 : length(individual(F(front).f(i)).p)

                 individual(individual(F(front).f(i)).p(j)).n = ...

                  individual(individual(F(front).f(i)).p(j)).n - 1;

                  if individual(individual(F(front).f(i)).p(j)).n == 0

                         x(individual(F(front).f(i)).p(j),M + V + 1) = ...

                       front + 1;

                   Q = [Q individual(F(front).f(i)).p(j)];

               end

           end

      end

   end

   front =  front + 1;

   F(front).f = Q;

end

[temp,index_of_fronts] = sort(x(:,M + V + 1));

for i = 1 : length(index_of_fronts)

    sorted_based_on_front(i,:) = x(index_of_fronts(i),:);

end

current_index = 0;

% Find the crowding distance for each individual in each front

for front = 1 : (length(F) - 1)

    objective = [];

    distance = 0;

    y = [];

    previous_index = current_index + 1;

    for i = 1 : length(F(front).f)

       y(i,:) = sorted_based_on_front(current_index + i,:);

    end

    current_index = current_index + i;

    % Sort each individual based on the objective

    sorted_based_on_objective = [];

    for i = 1 : M

       [sorted_based_on_objective, index_of_objectives] = ...

           sort(y(:,V + i));

       sorted_based_on_objective = [];

       for j = 1 : length(index_of_objectives)

           sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);

       end

       f_max = ...

           sorted_based_on_objective(length(index_of_objectives), V + i);

       f_min = sorted_based_on_objective(1, V + i);

       y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i)...

           = Inf;

       y(index_of_objectives(1),M + V + 1 + i) = Inf;

        for j = 2 : length(index_of_objectives) - 1

           next_obj  = sorted_based_on_objective(j + 1,V + i);

           previous_obj  = sorted_based_on_objective(j - 1,V + i);

           if (f_max - f_min == 0)

               y(index_of_objectives(j),M + V + 1 + i) = Inf;

           else

               y(index_of_objectives(j),M + V + 1 + i) = ...

                    (next_obj - previous_obj)/(f_max - f_min);

           end

        end

    end

    distance = [];

    distance(:,1) = zeros(length(F(front).f),1);

    for i = 1 : M

       distance(:,1) = distance(:,1) + y(:,M + V + 1 + i);

    end

    y(:,M + V + 2) = distance;

    y = y(:,1 : M + V + 2);

    z(previous_index:current_index,:) = y;

end

f = z();

Published with MATLAB® 7.0

Main Function

Main program to run the NSGA-II MOEA. Read thecorresponding documentation to learn more about multiobjective optimizationusing evolutionary algorithms. initialize_variables has two arguments; Firstbeing the population size and the second the problem number.
'1' corresponds toMOP1 and '2' corresponds to MOP2.

Contents

Initialize the variables

Sort the initialized population

Start the evolution process

Result

Visualize

Initialize the variables

Declare the variables and initialize their values pop -population gen - generations pro - problem number

pop = 200;

gen = 1;

pro = 1;

switch pro

    case 1

       % M is the number of objectives.

       M = 2;

       % V is the number of decision variables. In this case it is

       % difficult to visualize the decision variables space while the

       % objective space is just two dimensional.

       V = 6;

    case 2

       M = 3;

       V = 12;

end

% Initialize the population

chromosome = initialize_variables(pop,pro);

Sort the initialized population

Sort the population using non-domination-sort. This returnstwo columns for each individual which are the rank and the crowding distancecorresponding to their position in the front they belong.

chromosome = non_domination_sort_mod(chromosome,pro);

Start the evolution process

% The following are performed in each generation

% Select the parents

% Perfrom crossover and Mutation operator

% Perform Selection

for i = 1 : gen

    % Select the parents

    % Parents are selected for reproduction to generate offspring. The

    % original NSGA-II uses a binary tournament selection based on the

    % crowded-comparision operator. The arguments are

    % pool - size of the mating pool. It is common to have this to be half the

    %       population size.

    % tour - Tournament size. Original NSGA-II uses a binary tournament

    %       selection, but to see the effect of tournament size this is kept

    %       arbitary, to be choosen by the user.

    pool = round(pop/2);

    tour = 2;

    parent_chromosome = tournament_selection(chromosome,pool,tour);

    % Perfrom crossover and Mutation operator

    % The original NSGA-II algorithm uses Simulated Binary Crossover (SBX) and

    % Polynomial crossover. Crossover probability pc = 0.9 and mutation

    % probability is pm = 1/n, where n is the number of decision variables.

    % Both real-coded GA and binary-coded GA are implemented in the original

    % algorithm, while in this program only the real-coded GA is considered.

    % The distribution indeices for crossover and mutation operators as mu = 20

    % and mum = 20 respectively.

    mu = 20;

    mum = 20;

    offspring_chromosome = genetic_operator(parent_chromosome,pro,mu,mum);

    % Intermediate population

    % Intermediate population is the combined population of parents and

    % offsprings of the current generation. The population size is almost 1 and

    % half times the initial population.

    [main_pop,temp] = size(chromosome);

    [offspring_pop,temp] = size(offspring_chromosome);

    intermediate_chromosome(1:main_pop,:) = chromosome;

    intermediate_chromosome(main_pop + 1 : main_pop + offspring_pop,1 : M+V) = ...

       offspring_chromosome;

    % Non-domination-sort of intermediate population

    % The intermediate population is sorted again based on non-domination sort

    % before the replacement operator is performed on the intermediate

    % population.

    intermediate_chromosome = ...

       non_domination_sort_mod(intermediate_chromosome,pro);

    % Perform Selection

    % Once the intermediate population is sorted only the best solution is

    % selected based on it rank and crowding distance. Each front is filled in

    % ascending order until the addition of population size is reached. The

    % last front is included in the population based on the individuals with

    % least crowding distance

    chromosome = replace_chromosome(intermediate_chromosome,pro,pop);

    if ~mod(i,10)

       fprintf('%d\n',i);

    end

end

Result

Save the result in ASCII text format.

save solution.txt chromosome -ASCII

Visualize

The following is used to visualize the result for the givenproblem.

switch pro

    case 1

       plot(chromosome(:,V + 1),chromosome(:,V + 2),'*');

       title('MOP1 using NSGA-II');

       xlabel('f(x_1)');

       ylabel('f(x_2)');

    case 2

       plot3(chromosome(:,V + 1),chromosome(:,V + 2),chromosome(:,V + 3),'*');

       title('MOP2 using NSGA-II');

       xlabel('f(x_1)');

       ylabel('f(x_2)');

       zlabel('f(x_3)');

end

Published with MATLAB® 7.0

replace_chromosome(intermediate_chromosome,pro,pop)

This function replaces the chromosomes based on rank andcrowding distance. Initially until the population size is reached each front isadded one by one until addition of a complete front which results in exceedingthe population size. At this point the chromosomes
in that front is addedsubsequently to the population based on crowding distance.

[N,V] = size(intermediate_chromosome);

switch pro

       
d241
case 1

       M = 2;

       V = 6;

    case 2

       M = 3;

       V = 12;

end

% Get the index for the population sort based on the rank

[temp,index] = sort(intermediate_chromosome(:,M + V + 1));

% Now sort the individuals based on the index

for i = 1 : N

    sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);

end

% Find the maximum rank in the current population

max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the

% whole population is filled.

previous_index = 0;

for i = 1 : max_rank

    current_index = max(find(sorted_chromosome(:,M + V + 1) == i));

    if current_index > pop

       remaining = pop - previous_index;

       temp_pop = ...

           sorted_chromosome(previous_index + 1 : current_index, :);

       [temp_sort,temp_sort_index] = ...

           sort(temp_pop(:, M + V + 2),'descend');

       for j = 1 : remaining

           f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);

       end

       return;

    elseif current_index < pop

       f(previous_index + 1 : current_index, :) = ...

           sorted_chromosome(previous_index + 1 : current_index, :);

    else

       f(previous_index + 1 : current_index, :) = ...

           sorted_chromosome(previous_index + 1 : current_index, :);

       return;

    end

    previous_index = current_index;

end

Published with MATLAB® 7.0

function f =selection_individuals(chromosome,pool_size,tour_size)

% function selection_individuals(chromosome,pool_size,tour_size) is the

% selection policy for selecting the individuals for the mating pool. The

% selection is based on tournament selection. Argument 'chromosome' is the

% current generation population from which the individuals are selected to

% form a mating pool of size 'pool_size' after performing tournament

% selection, with size of the tournament being 'tour_size'. By varying the

% tournament size the selection pressure can be adjusted.

[pop,variables] = size(chromosome);

rank = variables - 1;

distance = variables;

for i = 1 : pool_size

    for j = 1 : tour_size

       candidate(j) = round(pop*rand(1));

       if candidate(j) == 0

           candidate(j) = 1;

       end

       if j > 1

           while ~isempty(find(candidate(1 : j - 1) == candidate(j)))

               candidate(j) = round(pop*rand(1));

               if candidate(j) == 0

                   candidate(j) = 1;

               end

           end

       end

    end

    for j = 1 : tour_size

       c_obj_rank(j) = chromosome(candidate(j),rank);

       c_obj_distance(j) = chromosome(candidate(j),distance);

    end

    min_candidate = ...

       find(c_obj_rank == min(c_obj_rank));

    if length(min_candidate) ~= 1

       max_candidate = ...

       find(c_obj_distance(min_candidate) == max(c_obj_distance(min_candidate)));

       if length(max_candidate) ~= 1

           max_candidate = max_candidate(1);

       end

       f(i,:) = chromosome(candidate(min_candidate(max_candidate)),:);

    else

       f(i,:) = chromosome(candidate(min_candidate(1)),:);

    end

end

Published with MATLAB® 7.0

 

function f = crowding_distance(x,problem)
% This function calculates the crowding distance
[N,M] = size(x);
switch problem
case 1
M = 2;
V = 6;
case 2
M = 3;
V = 12;
end

% Crowding distance for each front
for i = 1 : length(F(front).f)
y(i,:) = x(F(front).f(i),:);
end
for i = 1 : M
[sorted(i).individual,sorted(i).index] = sort(y(:,V + i));
distance(sorted(i).index(1)).individual = Inf;
distance(sorted(i).index(length(sorted(i).index))).individual = Inf;
end

[num,len] = size(y);
% Initialize all the distance of individuals as zero.
for i = 1 : M
for j = 2 : num - 1
distance(j).individual = 0;
end
objective(i).range = ...
sorted(i).individual(length(sorted(i).individual)) - ...
sorted(i).individual(1);
% Maximum and minimum objectives value for the ith objective
end
% Caluclate the crowding distance for front one.
for i = 1 : M
for j = 2 : num - 1
distance(j).individual = distance(j).individual + ...
(sorted(i).individual(j + 1) - sorted(i).individual(j - 1))/...
objective(i).range;
y(sorted(i).index(j),M + V + 2) = distance(j).individual;
end
end
</PRE

Published with MATLAB® 7.0
function f = genetic_operator(parent_chromosome,pro,mu,mum);
% This function is utilized to produce offsprings from parent chromosomes.
% The genetic operators corssover and mutation which are carried out with
% slight modifications from the original design. For more information read
% the document enclosed.

[N,M] = size(parent_chromosome);
switch pro
case 1
M = 2;
V = 6;
case 2
M = 3;
V = 12;
end
p = 1;
was_crossover = 0;
was_mutation = 0;
l_limit = 0;
u_limit = 1;
for i = 1 : N
if rand(1) < 0.9
child_1 = [];
child_2 = [];
parent_1 = round(N*rand(1));
if parent_1 < 1
parent_1 = 1;
end
parent_2 = round(N*rand(1));
if parent_2 < 1
parent_2 = 1;
end
while isequal(parent_chromosome(parent_1,:),parent_chromosome(parent_2,:))
parent_2 = round(N*rand(1));
if parent_2 < 1
parent_2 = 1;
end
end
parent_1 = parent_chromosome(parent_1,:);
parent_2 = parent_chromosome(parent_2,:);
for j = 1 : V
% SBX (Simulated Binary Crossover)
% Generate a random number
u(j) = rand(1);
if u(j) <= 0.5
bq(j) = (2*u(j))^(1/(mu+1));
else
bq(j) = (1/(2*(1 - u(j))))^(1/(mu+1));
end
child_1(j) = ...
0.5*(((1 + bq(j))*parent_1(j)) + (1 - bq(j))*parent_2(j));
child_2(j) = ...
0.5*(((1 - bq(j))*parent_1(j)) + (1 + bq(j))*parent_2(j));
if child_1(j) > u_limit
child_1(j) = u_limit;
elseif child_1(j) < l_limit
child_1(j) = l_limit;
end
if child_2(j) > u_limit
child_2(j) = u_limit;
elseif child_2(j) < l_limit
child_2(j) = l_limit;
end
end
child_1(:,V + 1: M + V) = evaluate_objective(child_1,pro);
child_2(:,V + 1: M + V) = evaluate_objective(child_2,pro);
was_crossover = 1;
was_mutation = 0;
else
parent_3 = round(N*rand(1));
if parent_3 < 1
parent_3 = 1;
end
% Make sure that the mutation does not result in variables out of
% the search space. For both the MOP's the range for decision space
% is [0,1]. In case different variables have different decision
% space each variable can be assigned a range.
child_3 = parent_chromosome(parent_3,:);
for j = 1 : V
r(j) = rand(1);
if r(j) < 0.5
delta(j) = (2*r(j))^(1/(mum+1)) - 1;
else
delta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1));
end
child_3(j) = child_3(j) + delta(j);
if child_3(j) > u_limit
child_3(j) = u_limit;
elseif child_3(j) < l_limit
child_3(j) = l_limit;
end
end
child_3(:,V + 1: M + V) = evaluate_objective(child_3,pro);
was_mutation = 1;
was_crossover = 0;
end
if was_crossover
child(p,:) = child_1;
child(p+1,:) = child_2;
was_cossover = 0;
p = p + 2;
elseif was_mutation
child(p,:) = child_3(1,1 : M + V);
was_mutation = 0;
p = p + 1;
end
end
f = child;

Published with MATLAB® 7.0
function f = initialize_variables(N,problem)
% function f = initialize_variables(N,problem)
% N - Population size
% problem - takes integer values 1 and 2 where,
%           '1' for MOP1
%           '2' for MOP2
%
% This function initializes the population with N individuals and each
% individual having M decision variables based on the selected problem.
% M = 6 for problem MOP1 and M = 12 for problem MOP2. The objective space
% for MOP1 is 2 dimensional while for MOP2 is 3 dimensional.

% Both the MOP's has 0 to 1 as its range for all the decision variables.
min = 0;
max = 1;
switch problem
case 1
M = 6;
K = 8;
case 2
M = 12;
K = 15;
end
for i = 1 : N
% Initialize the decision variables
for j = 1 : M
f(i,j) = rand(1); % i.e f(i,j) = min + (max - min)*rand(1);
end
% Evaluate the objective function
f(i,M + 1: K) = evaluate_objective(f(i,:),problem);
end
</PRE

Published with MATLAB® 7.0

Non-Donimation Sort
This function sort the current popultion based on non-domination. All the individuals in the first front are given a rank of 1, the second front individuals are assigned rank 2 and so on. After assigning the rank the crowding in each front is calculated.
[N,M] = size(x);
switch problem
case 1
M = 2;
V = 6;
case 2
M = 3;
V = 12;
end
front = 1;

% There is nothing to this assignment, used only to manipulate easily in
% MATLAB.
F(front).f = [];
individual = [];
for i = 1 : N
% Number of individuals that dominate this individual
individual(i).n = 0;
% Individuals which this individual dominate
individual(i).p = [];
for j = 1 : N
dom_less = 0;
dom_equal = 0;
dom_more = 0;
for k = 1 : M
if (x(i,V + k) < x(j,V + k))
dom_less = dom_less + 1;
elseif (x(i,V + k) == x(j,V + k))
dom_equal = dom_equal + 1;
else
dom_more = dom_more + 1;
end
end
if dom_less == 0 & dom_equal ~= M
individual(i).n = individual(i).n + 1;
elseif dom_more == 0 & dom_equal ~= M
individual(i).p = [individual(i).p j];
end
end
if individual(i).n == 0
x(i,M + V + 1) = 1;
F(front).f = [F(front).f i];
end
end
% Find the subsequent fronts
while ~isempty(F(front).f)
Q = [];
for i = 1 : length(F(front).f)
if ~isempty(individual(F(front).f(i)).p)
for j = 1 : length(individual(F(front).f(i)).p)
individual(individual(F(front).f(i)).p(j)).n = ...
individual(individual(F(front).f(i)).p(j)).n - 1;
if individual(individual(F(front).f(i)).p(j)).n == 0
x(individual(F(front).f(i)).p(j),M + V + 1) = ...
front + 1;
Q = [Q individual(F(front).f(i)).p(j)];
end
end
end
end
front =  front + 1;
F(front).f = Q;
end
[temp,index_of_fronts] = sort(x(:,M + V + 1));
for i = 1 : length(index_of_fronts)
sorted_based_on_front(i,:) = x(index_of_fronts(i),:);
end
current_index = 0;
% Find the crowding distance for each individual in each front
for front = 1 : (length(F) - 1)
objective = [];
distance = 0;
y = [];
previous_index = current_index + 1;
for i = 1 : length(F(front).f)
y(i,:) = sorted_based_on_front(current_index + i,:);
end
current_index = current_index + i;
% Sort each individual based on the objective
sorted_based_on_objective = [];
for i = 1 : M
[sorted_based_on_objective, index_of_objectives] = ...
sort(y(:,V + i));
sorted_based_on_objective = [];
for j = 1 : length(index_of_objectives)
sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);
end
f_max = ...
sorted_based_on_objective(length(index_of_objectives), V + i);
f_min = sorted_based_on_objective(1, V + i);
y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i)...
= Inf;
y(index_of_objectives(1),M + V + 1 + i) = Inf;
for j = 2 : length(index_of_objectives) - 1
next_obj  = sorted_based_on_objective(j + 1,V + i);
previous_obj  = sorted_based_on_objective(j - 1,V + i);
if (f_max - f_min == 0)
y(index_of_objectives(j),M + V + 1 + i) = Inf;
else
y(index_of_objectives(j),M + V + 1 + i) = ...
(next_obj - previous_obj)/(f_max - f_min);
end
end
end
distance = [];
distance(:,1) = zeros(length(F(front).f),1);
for i = 1 : M
distance(:,1) = distance(:,1) + y(:,M + V + 1 + i);
end
y(:,M + V + 2) = distance;
y = y(:,1 : M + V + 2);
z(previous_index:current_index,:) = y;
end
f = z();

Published with MATLAB® 7.0
Main Function
Main program to run the NSGA-II MOEA. Read the corresponding documentation to learn more about multiobjective optimization using evolutionary algorithms. initialize_variables has two arguments; First being the population size and the second the problem number. '1' corresponds to MOP1 and '2' corresponds to MOP2.
Contents
•	Initialize the variables
•	Sort the initialized population
•	Start the evolution process
•	Result
•	Visualize
Initialize the variables
Declare the variables and initialize their values pop - population gen - generations pro - problem number
pop = 200;
gen = 1;
pro = 1;

switch pro
case 1
% M is the number of objectives.
M = 2;
% V is the number of decision variables. In this case it is
% difficult to visualize the decision variables space while the
% objective space is just two dimensional.
V = 6;
case 2
M = 3;
V = 12;
end

% Initialize the population
chromosome = initialize_variables(pop,pro);
Sort the initialized population
Sort the population using non-domination-sort. This returns two columns for each individual which are the rank and the crowding distance corresponding to their position in the front they belong.
chromosome = non_domination_sort_mod(chromosome,pro);
Start the evolution process
% The following are performed in each generation
% Select the parents
% Perfrom crossover and Mutation operator
% Perform Selection

for i = 1 : gen
% Select the parents
% Parents are selected for reproduction to generate offspring. The
% original NSGA-II uses a binary tournament selection based on the
% crowded-comparision operator. The arguments are
% pool - size of the mating pool. It is common to have this to be half the
%        population size.
% tour - Tournament size. Original NSGA-II uses a binary tournament
%        selection, but to see the effect of tournament size this is kept
%        arbitary, to be choosen by the user.
pool = round(pop/2);
tour = 2;
parent_chromosome = tournament_selection(chromosome,pool,tour);

% Perfrom crossover and Mutation operator
% The original NSGA-II algorithm uses Simulated Binary Crossover (SBX) and
% Polynomial crossover. Crossover probability pc = 0.9 and mutation
% probability is pm = 1/n, where n is the number of decision variables.
% Both real-coded GA and binary-coded GA are implemented in the original
% algorithm, while in this program only the real-coded GA is considered.
% The distribution indeices for crossover and mutation operators as mu = 20
% and mum = 20 respectively.
mu = 20;
mum = 20;
offspring_chromosome = genetic_operator(parent_chromosome,pro,mu,mum);

% Intermediate population
% Intermediate population is the combined population of parents and
% offsprings of the current generation. The population size is almost 1 and
% half times the initial population.
[main_pop,temp] = size(chromosome);
[offspring_pop,temp] = size(offspring_chromosome);
intermediate_chromosome(1:main_pop,:) = chromosome;
intermediate_chromosome(main_pop + 1 : main_pop + offspring_pop,1 : M+V) = ...
offspring_chromosome;

% Non-domination-sort of intermediate population
% The intermediate population is sorted again based on non-domination sort
% before the replacement operator is performed on the intermediate
% population.
intermediate_chromosome = ...
non_domination_sort_mod(intermediate_chromosome,pro);
% Perform Selection
% Once the intermediate population is sorted only the best solution is
% selected based on it rank and crowding distance. Each front is filled in
% ascending order until the addition of population size is reached. The
% last front is included in the population based on the individuals with
% least crowding distance
chromosome = replace_chromosome(intermediate_chromosome,pro,pop);
if ~mod(i,10)
fprintf('%d\n',i);
end
end
Result
Save the result in ASCII text format.
save solution.txt chromosome -ASCII
Visualize
The following is used to visualize the result for the given problem.
switch pro
case 1
plot(chromosome(:,V + 1),chromosome(:,V + 2),'*');
title('MOP1 using NSGA-II');
xlabel('f(x_1)');
ylabel('f(x_2)');
case 2
plot3(chromosome(:,V + 1),chromosome(:,V + 2),chromosome(:,V + 3),'*');
title('MOP2 using NSGA-II');
xlabel('f(x_1)');
ylabel('f(x_2)');
zlabel('f(x_3)');
end

Published with MATLAB® 7.0
replace_chromosome(intermediate_chromosome,pro,pop)
This function replaces the chromosomes based on rank and crowding distance. Initially until the population size is reached each front is added one by one until addition of a complete front which results in exceeding the population size. At this point the chromosomes in that front is added subsequently to the population based on crowding distance.
[N,V] = size(intermediate_chromosome);
switch pro
case 1
M = 2;
V = 6;
case 2
M = 3;
V = 12;
end

% Get the index for the population sort based on the rank
[temp,index] = sort(intermediate_chromosome(:,M + V + 1));

% Now sort the individuals based on the index
for i = 1 : N
sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end

% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
current_index = max(find(sorted_chromosome(:,M + V + 1) == i));
if current_index > pop
remaining = pop - previous_index;
temp_pop = ...
sorted_chromosome(previous_index + 1 : current_index, :);
[temp_sort,temp_sort_index] = ...
sort(temp_pop(:, M + V + 2),'descend');
for j = 1 : remaining
f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
end
return;
elseif current_index < pop
f(previous_index + 1 : current_index, :) = ...
sorted_chromosome(previous_index + 1 : current_index, :);
else
f(previous_index + 1 : current_index, :) = ...
sorted_chromosome(previous_index + 1 : current_index, :);
return;
end
previous_index = current_index;
end

Published with MATLAB® 7.0
function f = selection_individuals(chromosome,pool_size,tour_size)
% function selection_individuals(chromosome,pool_size,tour_size) is the
% selection policy for selecting the individuals for the mating pool. The
% selection is based on tournament selection. Argument 'chromosome' is the
% current generation population from which the individuals are selected to
% form a mating pool of size 'pool_size' after performing tournament
% selection, with size of the tournament being 'tour_size'. By varying the
% tournament size the selection pressure can be adjusted.

[pop,variables] = size(chromosome);
rank = variables - 1;
distance = variables;

for i = 1 : pool_size
for j = 1 : tour_size
candidate(j) = round(pop*rand(1));
if candidate(j) == 0
candidate(j) = 1;
end
if j > 1
while ~isempty(find(candidate(1 : j - 1) == candidate(j)))
candidate(j) = round(pop*rand(1));
if candidate(j) == 0
candidate(j) = 1;
end
end
end
end
for j = 1 : tour_size
c_obj_rank(j) = chromosome(candidate(j),rank);
c_obj_distance(j) = chromosome(candidate(j),distance);
end
min_candidate = ...
find(c_obj_rank == min(c_obj_rank));
if length(min_candidate) ~= 1
max_candidate = ...
find(c_obj_distance(min_candidate) == max(c_obj_distance(min_candidate)));
if length(max_candidate) ~= 1
max_candidate = max_candidate(1);
end
f(i,:) = chromosome(candidate(min_candidate(max_candidate)),:);
else
f(i,:) = chromosome(candidate(min_candidate(1)),:);
end
end

Published with MATLAB® 7.0