NSGA_II学习案例1

此文章来自网上，转来学习一下。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% function nsga_2(pro)
%% 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.
%inp_para_definition=input_parameters_definition;
%% Initialize the variables
% Declare the variables and initialize their values
% pop - population
% gen - generations
% pro - problem number
%clear;clc;tic;
pop = 100; % 每一代的种群数
gen = 100; % 总共的代数
pro = 2; % 问题选择1或者2，见switch
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;
case 3 % case 1和case 2 用来对整个算法进行常规验证，作为调试之用；case 3 为本工程所需；
M = 2; %(output parameters 个数)
V = 8; %(input parameters 个数)
K = 10;
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. 真是牛 X了。
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;
% 针对对象是上一步产生的新的个体parent_chromosome
%对parent_chromosome 每次操作以较大的概率进行交叉（产生两个新的候选人），或者较小的概率变异（一个新的候选人）操作，这样
%就会产生较多的新个体
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;
%intermediate_chromosome=inter_chromo(chromosome,offspring_chromosome,pro);


% 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),’y+’);
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
%disp(‘run time is:’)
%toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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; % k=决策变量(M=6)+目标变量(K-M=2)=8
case 2
M = 12;
K = 15;
case 3 % case 1和case 2 用来对整个算法进行常规验证，作为调试之用；case 3 为本工程所需；
M = 8; %(input parameters 个数)
K = 10;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = evaluate_objective(x,problem)
% Function to evaluate the objective functions for the given input vector
% x. x has the decision variables
switch problem
case 1
f = [];
%% Objective function one
f(1) = 1 - exp(-4*x(1))*(sin(6*pi*x(1)))^6;
sum = 0;
for i = 2 : 6
sum = sum + x(i)/4;
end
%% Intermediate function
g_x = 1 + 9*(sum)^(0.25);
%% Objective function one
f(2) = g_x*(1 - ((f(1))/(g_x))^2);
case 2
f = [];
%% Intermediate function
g_x = 0;
for i = 3 : 12
g_x = g_x + (x(i) - 0.5)^2;
end
%% Objective function one
f(1) = (1 + g_x)*cos(0.5*pi*x(1))*cos(0.5*pi*x(2));
%% Objective function two
f(2) = (1 + g_x)*cos(0.5*pi*x(1))*sin(0.5*pi*x(2));
%% Objective function three
f(3) = (1 + g_x)*sin(0.5*pi*x(1));
case 3
f = [];
%% Objective function one
f(1) = 0;
%% Objective function one
f(2) = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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.
function f = non_domination_sort_mod(x,problem)
[N,M] = size(x);
switch problem
case 1
M = 2;
V = 6;
case 2
M = 3;
V = 12;
case 3 % case 1和case 2 用来对整个算法进行常规验证，作为调试之用；case 3 为本工程所需；
M = 2; %(output parameters 个数)
V = 8; %(input parameters 个数)
K = 10;
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 支配i的解的个数
individual(i).n = 0;
% Individuals which this individual dominate 被i支配的解
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 % 举个例子,其中i=a,b; j= c,d; 如果i中的a,b全部大于
individual(i).n = individual(i).n + 1; % 或者部分大于j中的c,d(但没有小于的情况),则称为i优于j,
elseif dom_more == 0 & dom_equal ~= M % 当i优于j的时候,则此时把individual(i)_n加1
individual(i).p = [individual(i).p j]; %如果i中的a,b全部小于或者部分小于j中的c,d(但没有大于的情况),则
end %则称为j优于i, 则把此时的j放入individual(i)_P中;
end %总之，就是说两个目标变量必须全部大于或者全部小于才能对individual有效。
if individual(i).n == 0 % 如果没有劣于i的话,即F(front).f 存放的都是差的解.
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) %i表示最优解
if ~isempty(individual(F(front).f(i)).p)
for j = 1 : length(individual(F(front).f(i)).p) % 被前端解i所支配的解得集合，这个集合由j构成，
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
% function sort : sort(x_sequence)===>>return a increasing orer data sequence
[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),:); % 对解（染色体）进行排序，按照front分层
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; % 最大值放lnf
y(index_of_objectives(1),M + V + 1 + i) = Inf; % 最小值放lnf
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = replace_chromosome(intermediate_chromosome,pro,pop)
%% 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;
case 3 % case 1和case 2 用来对整个算法进行常规验证，作为调试之用；case 3 为本工程所需；
M = 2; %(output parameters 个数)
V = 8; %(input parameters 个数)
K = 10;
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_pop_1=temp_pop*(-1);
[temp_sort_1,temp_sort_index] = …
sort(temp_pop_1(:, M + V + 2)); %在编译的时候出现的问题找到了，主要是对sort(a,’descend’)并不支持。
temp_sort=temp_sort_1*(-1);
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = tournament_selection(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 % 这个loop目的是为了随机找出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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
%input_parameters=inp_para_definition; %定义了一个输入变量的结构,
[N,M] = size(parent_chromosome); % 注意这里的M其实是个垃圾，后面语句会对其重新赋值利用
switch pro
case 1
M = 2;
V = 6;
case 2
M = 3;
V = 12;
case 3 % case 1和case 2 用来对整个算法进行常规验证，作为调试之用；case 3 为本工程所需；
M = 2; %(output parameters 个数)
V = 8; %(input parameters 个数)
K = 10;
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%