目录
一、NSGA-II 算法流程图
二、部分函数详细注释
1、主函数(nsga_2_optimization)
2、初始化代码 (initialize_variables)
3、快速非支配排序和拥挤度计算(non_domination_sort_mod)
4、生成新的种群、精英策略(replace_chromosome)
5、目标函数(evaluate_objective)
一、NSGA-II 算法流程图
二、部分函数详细注释
1、主函数(nsga_2_optimization)
function nsga_2_optimization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pop = 200; %种群数量
gen = 500; %迭代次数
M = 2; %目标函数数量
V = 30; %维度(决策变量的个数) 决策变量就是解的个数
min_range = zeros(1, V); %下界 生成1*30的个体向量 全为0
max_range = ones(1,V); %上界 生成1*30的个体向量 全为1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
chromosome = initialize_variables(pop, M, V, min_range, max_range);%初始化种群
chromosome = non_domination_sort_mod(chromosome, M, V);%对初始化种群进行非支配快速排序和拥挤度计算
for i = 1 : gen
pool = round(pop/2);%round() 四舍五入取整 交配池大小
tour = 2;%竞标赛 参赛选手个数
parent_chromosome = tournament_selection(chromosome, pool, tour);%竞标赛选择适合繁殖的父代
mu = 20;%交叉和变异算法的分布指数
mum = 20;
%% parent_chromosome 竞标赛选择的适合繁殖的父代 M 目标函数数量 V维度(决策变量的个数) mu = 20;%交叉和变异算法的分布指数 mum = 20;
%%min_range = zeros(1, V); %下界 生成1*30的个体向量 全为0
%%max_range = ones(1,V); %上界 生成1*30的个体向量 全为1 这个在这里用来约束解(Xi)的范围
%%offspring_chromosome不一定生成了500个后代
offspring_chromosome = genetic_operator(parent_chromosome,M, V, mu, mum, min_range, max_range);%进行交叉变异产生子代 该代码中使用模拟二进制交叉和多项式变异 采用实数编码
[main_pop,~] = size(chromosome);%父代种群的大小
[offspring_pop,~] = size(offspring_chromosome);%子代种群的大小
clear temp
intermediate_chromosome(1:main_pop,:) = chromosome;
intermediate_chromosome(main_pop + 1 : main_pop + offspring_pop,1 : M+V) = offspring_chromosome;%合并父代种群和子代种群
intermediate_chromosome = non_domination_sort_mod(intermediate_chromosome, M, V);%对新的种群进行快速非支配排序
%%精英选择 从子代和父代中选出Pop个
chromosome = replace_chromosome(intermediate_chromosome, M, V, pop);%选择合并种群中前N个优先的个体组成新种群
%%每计算100代清空下控制台
if ~mod(i,100)
clc;
fprintf('%d generations completed\n',i);
end
end
if M == 2
plot(chromosome(:,V + 1),chromosome(:,V + 2),'*');
xlabel('f_1'); ylabel('f_2');
title('Pareto Optimal Front');
elseif M == 3
plot3(chromosome(:,V + 1),chromosome(:,V + 2),chromosome(:,V + 3),'*');
xlabel('f_1'); ylabel('f_2'); zlabel('f_3');
title('Pareto Optimal Surface');
end
2、初始化代码 (initialize_variables)
function f = initialize_variables(N, M, V, min_range, max_range)%f是一个由种群个体组成的矩阵
%M目标函数数量
%V维度(决策变量的个数)
%N 种群数量
min = min_range; %下界 生成1*30的个体向量 全为0
max = max_range; %上界 生成1*30的个体向量 全为1
K = M + V;%%K是数组的总元素个数。为了便于计算,决策变量和目标函数串在一起形成一个数组。
%对于交叉和变异,利用目标变量对决策变量进行选择
for i = 1 : N
for j = 1 : V
f(i,j) = min(j) + (max(j) - min(j))*rand(1);%f(i j)表示的是种群中第i个个体中的第j个决策变量, %%-ppppppppppp
%这行代码为每个个体的所有决策变量在约束条件内随机取值
end
f(i,V + 1: K) = evaluate_objective(f(i,:), M, V); % M是目标函数数量 V是决策变量个数 f(:,1)就是取f 矩阵的第1列。
%为了简化计算将对应的目标函数值储存在染色体的V + 1 到 K的位置。
%%k+1-V 存的是目标函数的值
end
3、快速非支配排序和拥挤度计算(non_domination_sort_mod)
%% 对初始种群开始排序 快速非支配排序
% 使用非支配排序对种群进行排序。该函数返回每个个体对应的排序值和拥挤距离,是一个两列的矩阵。
% 并将排序值和拥挤距离添加到染色体矩阵中
%M + V + 1 M + V + 2 存的是分层等级
%x 进来时是200行(种群个数)*32列(30列x1-x30 2列 f1 f2)
%x出去时33列存的是当前个体的层级 1为最高
function f = non_domination_sort_mod(x, M, V)
%%chromosome = non_domination_sort_mod(chromosome, M, V); x/chromosome是解
%%其中包括了30个解和 目标函数的值 共32个 M 目标函数数量 V决策变量个数
[N, ~] = size(x);% N为矩阵x的行数,也是种群的数量
clear m
front = 1;%front记录了当前正在筛选那一层级的个体
F(front).f = []; %f是一个数组 用于存放当前front层级的个体
individual = [];%%存放个体i的支配个体的信息
for i = 1 : N
individual(i).n = 0;%n是个体i被支配的个体数量
individual(i).p = [];%p是被个体i支配的个体集合
for j = 1 : N
dom_less = 0;
dom_equal = 0;
dom_more = 0;
for k = 1 : M %判断个体i和个体j的支配关系
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受j支配,相应的n加1
individual(i).n = individual(i).n + 1;
elseif dom_more == 0 && dom_equal ~= M % 说明i支配j,把j加入i的支配合集中
individual(i).p = [individual(i).p j];
end
end
%找出最高等级的所有个体
if individual(i).n == 0 %个体i非支配等级排序最高,属于当前最优解集,相应的染色体中携带代表排序数的信息
x(i,M + V + 1) = 1;%1代表最高等级 改个体i的层级为1
F(front).f = [F(front).f i];%等级为1的非支配解集 f是个矩阵 在F(front).f 矩阵后面加上 i 赋值给F(front).f
end
end
%上面的代码是为了找出等级最高的非支配解集
%下面的代码是为了给其他个体进行分级
while ~isempty(F(front).f)
Q = []; %存放下一个front集合
for i = 1 : length(F(front).f)%循环当前支配解集中的个体
if ~isempty(individual(F(front).f(i)).p)%个体i有自己所支配的解集
for j = 1 : length(individual(F(front).f(i)).p)%循环个体i所支配解集中的个体
%%individual(F(front).f(i)).p(j) 代表front层个体所支配的一个个体
individual(individual(F(front).f(i)).p(j)).n = ...%...表示的是与下一行代码是相连的, 这里表示个体j的被支配个数减1
individual(individual(F(front).f(i)).p(j)).n - 1; %因为层级最高为1(层级为1 即这时n为0)的在上面已经筛选完 这个循环里面的n最少为1 都被支配
%代表去掉front层的个体后,front层个体所支配的一个个体的被支配的个数为0时,这时它是这时候的优解之一
if individual(individual(F(front).f(i)).p(j)).n == 0% 如果q是非支配解集,则放入集合Q中
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
%x(:,M + V + 1)就是取x矩阵的第M + V + 1列 temp在这里没用 为了就是得到index_of_fronts
%这里想让这个矩阵按照第M+V+1也就层级进行排序,因为matlab没有直接排序的方法,所以采取这样的方法 x本身的顺序没有变
%sorted_based_on_front存储了排序后的矩阵
[temp,index_of_fronts] = sort(x(:,M + V + 1));%对个体的代表排序等级的列向量进行升序排序 index_of_fronts表示排序后的值对应原来的未排序的矩阵在当前列的行下标
%x(:,M + V + 1) 就是取这个矩阵的M + V + 1列
for i = 1 : length(index_of_fronts)
sorted_based_on_front(i,:) = x(index_of_fronts(i),:);%sorted_based_on_front中存放的是x矩阵按照排序等级升序排序后的矩阵 index_of_fronts(i)放到第i行
end
current_index = 0;
%% Crowding distance 计算每个个体的拥挤度
for front = 1 : (length(F) - 1)%这里减1是因为代码55行这里,F的最后一个元素为空,这样才能跳出循环。所以一共有length-1个排序等级
distance = 0;
y = [];
previous_index = current_index + 1; %% current_index记录的是每一层级的第一个个体在sorted_based_on_front的下标-1的值,用于每次
%循环,找到该层级的第一个个体
for i = 1 : length(F(front).f)
y(i,:) = sorted_based_on_front(current_index + i,:);%y中存放的是排序等级为front的集合矩阵
end
current_index = current_index + i;%current_index =i
sorted_based_on_objective = [];%存放基于拥挤距离排序的矩阵
for i = 1 : M
[sorted_based_on_objective, index_of_objectives] = ...
sort(y(:,V + i));%按照目标函数值排序 先按目标f1的值排序 再按f2排序
sorted_based_on_objective = [];
for j = 1 : length(index_of_objectives)
sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);% sorted_based_on_objective存放按照目标函数值排序后的y矩阵 在这里的y已经对层级进行过排序
end
f_max = ...
sorted_based_on_objective(length(index_of_objectives), V + i);%fmax为目标函数最大值 fmin为目标函数最小值
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) %%如果这一层级里只有一个个体 则这一层的个体fi的距离为无穷
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); %%这里是归一化处理 用最简单的标准化方法 让间距的值小一点 y的M + V + 1 + 1记录f1的拥挤度 M + V + 1 + 2记录
end
end
end
distance = [];
distance(:,1) = zeros(length(F(front).f),1);%%初始化 这个列向量用于保存当前front层的每一个体的拥挤度
for i = 1 : M
distance(:,1) = distance(:,1) + y(:,M + V + 1 + i); %%将y的f1的距离和f2的距离相加作为该个体的拥挤度
end
y(:,M + V + 2) = distance;
y = y(:,1 : M + V + 2); %%应该是可省的
z(previous_index:current_index,:) = y; %%在这里 previous_index记录该层级中第一个个体在矩阵中的行下标 current_index记录该层级中最后一个个体在矩阵中的行下标
end
f = z();%得到的是已经包含等级和拥挤度的种群矩阵 并且已经按等级和拥挤距离排序
4、生成新的种群、精英策略(replace_chromosome)
function f = replace_chromosome(intermediate_chromosome, M, V,pop)%精英选择策略
%%replace_chromosome(intermediate_chromosome, M, V, pop)
%%intermediate_chromosom 子代和父代两代种群放在一起进行非支配排序后的矩阵
%%生成新的种群(精英策略)
[N, m] = size(intermediate_chromosome);
%%获得按层级进行排序的种群索引
[temp,index] = sort(intermediate_chromosome(:,M + V + 1)); %%Matlab中没有办法直接让矩阵按照某一列直接排序所以只养血
clear temp m
%%获得按层级排序后的种群 前面非支配选择不是已经排过序了?
for i = 1 : N
sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end
%获得这个种群最大的层级
max_rank = max(intermediate_chromosome(:,M + V + 1));
%%用于保存当前新的种群中已经筛选到的个体数量
previous_index = 0;
for i = 1 : max_rank
%%找到当前层级的所有个体
current_index = max(find(sorted_chromosome(:,M + V + 1) == i));%%这样写避免了遍历寻找当前层级为i的个体 直接得到当前层级个体最大的索引,因为已经排过序
%%所以根据current_index就可以得到当前层级的所有个体
%%current_inde的位置就代表了当前已经筛选出了多少个个体了 因为前面排过序
%当该层级的个体大于所需的种群大小时
if current_index > pop
%%总共需要筛选出pop个个体,remaining保存还需筛选出多少个体新一代种群才达到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');
%%获得该层级中拥挤距离大的remaining个 个体
for j = 1 : remaining
f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
end
return;
elseif current_index < pop
%%小于Pop 说明新种群中个体数量<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
5、目标函数(evaluate_objective)
%%目标函数 ZDT1是MOP中常用的测试函数
%%这里有两个目标函数 f1 f2 这里就是求f1 f2值得公式
function f = evaluate_objective(x, M, V)%%计算每个个体的M个目标函数值
%%决策变量就是解的个数
%%f(i,V + 1: K) = evaluate_objective(f(i,:), M, V); % M是目标函数数量 V是决策变量个数 f(:,1)就是取f 矩阵的第1列。
f = [];
f(1) = x(1);
g = 1;
sum = 0;
for i = 2:V
sum = sum + x(i);
end
sum = 9*(sum / (V-1));
g = g + sum;
f(2) = g * (1 - sqrt(x(1) / g));
end