NSGA-II改进之实现代码（MATLAB）
引言
NSGA-II的改进内容：种群初始化、父代选择、种群交叉、种群变异。
具体的改进内容可查看：
改进前与改进后的优化结果的对比分析内容，以及对应的分析代码可查看下：
NSGA-II改进前后数据分析代码（matlab）
主要内容
此部分框架与多目标优化NSGA-II的实现（MATLAB完整代码）一致
区别在于：Mynsga2_test，Mynsga2_main，init_pop，myga，select_parent，这些内容上
代码模块
- my_nsga2_test：测试函数，用于保存测试数据
- my_nsga2_main：主函数，，用于运行NSGA2算法的框架
- get_variable_bounds：获取种群范围
- init_pop：种群初始化，采用佳点集方式生成种群
- sort_pop：种群排序
- select_parent：选择父代，采用基于线性排名的父代选择方式
- myga：进行遗传算法，杂交变异，采用非均匀变异方式和引导式交叉方式
- combined_pop：子代和原始种群进行合并
- select_pop：选择新一代种群
- calculate_gd：计算GD
- calculate_sp：计算SP
- calculate_pop：计算种群
其他内容
plotPareto：画出已知前言数据，用于跟测试得到的前言的可视化对比
如果需要获取已知的ZDT1、ZDT2、ZDT3、ZDT4、ZDT6的前言数据，通过以下链接获取：
ZDT前沿数据.zip
运行注意事项
在my_nsga2_test中设置pop_size，iterations。以及测试次数test。nsga2_test可以保存每一个测试函数，每一测试中每一代的种群数据以及GD和SP的数据。主要是为了方便用获取得到的数据进行分析。
如果只是需要查看效果，运行my_nsga2_main函数，注意pop_size，iterations的设置。
代码内容
my_nsga2_test
clc
clear
test = 1;%测试1次
% popsize，iterations在my_nsga2_main中也设置了，保持一致，可以设置为全局变量，未修改
popsize = 100;
iterations = 50;
for x = 1:6
switch x
case 1
[~,dim] = get_variable_bounds(x);
%dim = 10;
%设置保存参数
testGD = zeros(test,iterations);
testSP = zeros(test,iterations);
testPop = zeros(test,iterations,popsize,dim+4);
MyNSGA2_zdt1 = [];
disp('正在测试zdt1')
for i = 1:test
disp(['第' num2str(i) '次测试']);
[pop,GD,SP] = my_nsga2_main(x);
testGD(i,:) = GD;
testSP(i,:) = SP;
testPop(i,:,:,:) = pop;
end
MyNSGA2_zdt1.testGD = testGD;
MyNSGA2_zdt1.testSP = testSP;
MyNSGA2_zdt1.testPop = testPop;
save('MyNSGA2_zdt1.mat','MyNSGA2_zdt1')
case 2
[~,dim] = get_variable_bounds(x);
%dim = 10;
%设置保存参数
testGD = zeros(test,iterations);
testSP = zeros(test,iterations);
testPop = zeros(test,iterations,popsize,dim+4);
MyNSGA2_zdt2 = [];
disp('正在测试zdt2')
for i = 1:test
disp(['第' num2str(i) '次测试']);
[pop,GD,SP] = my_nsga2_main(x);
testGD(i,:) = GD;
testSP(i,:) = SP;
testPop(i,:,:,:) = pop;
end
MyNSGA2_zdt2.testGD = testGD;
MyNSGA2_zdt2.testSP = testSP;
MyNSGA2_zdt2.testPop = testPop;
save('MyNSGA2_zdt2.mat','MyNSGA2_zdt2')
case 3
[~,dim] = get_variable_bounds(x);
%dim = 10;
%设置保存参数
testGD = zeros(test,iterations);
testSP = zeros(test,iterations);
testPop = zeros(test,iterations,popsize,dim+4);
MyNSGA2_zdt3 = [];
disp('正在测试zdt3')
for i = 1:test
disp(['第' num2str(i) '次测试']);
[pop,GD,SP] = my_nsga2_main(x);
testGD(i,:) = GD;
testSP(i,:) = SP;
testPop(i,:,:,:) = pop;
end
MyNSGA2_zdt3.testGD = testGD;
MyNSGA2_zdt3.testSP = testSP;
MyNSGA2_zdt3.testPop = testPop;
save('MyNSGA2_zdt3.mat','MyNSGA2_zdt3')
case 4
[~,dim] = get_variable_bounds(x);
%dim = 10;
%设置保存参数
testGD = zeros(test,iterations);
testSP = zeros(test,iterations);
testPop = zeros(test,iterations,popsize,dim+4);
MyNSGA2_zdt4 = [];
disp('正在测试zdt4')
for i = 1:test
disp(['第' num2str(i) '次测试']);
[pop,GD,SP] = my_nsga2_main(x);
testGD(i,:) = GD;
testSP(i,:) = SP;
testPop(i,:,:,:) = pop;
end
MyNSGA2_zdt4.testGD = testGD;
MyNSGA2_zdt4.testSP = testSP;
MyNSGA2_zdt4.testPop = testPop;
save('MyNSGA2_zdt4.mat','MyNSGA2_zdt4')
case 5
[~,dim] = get_variable_bounds(x);
%dim = 10;
%设置保存参数
testGD = zeros(test,iterations);
testSP = zeros(test,iterations);
testPop = zeros(test,iterations,popsize,dim+4);
MyNSGA2_zdt6 = [];
disp('正在测试zdt6')
for i = 1:test
disp(['第' num2str(i) '次测试']);
[pop,GD,SP] = my_nsga2_main(x);
testGD(i,:) = GD;
testSP(i,:) = SP;
testPop(i,:,:,:) = pop;
end
MyNSGA2_zdt6.testGD = testGD;
MyNSGA2_zdt6.testSP = testSP;
MyNSGA2_zdt6.testPop = testPop;
save('MyNSGA2_zdt6.mat','MyNSGA2_zdt6')
end
end
my_nsga2_main
function [allpop,GD,SP] = my_nsga2_main(x)
% 测试主函数 x,问题编号
% 输出种群，GD和SP
% 参数设置 ,
% global pop_size
% global iterations;%迭代次数
target = 2;
pop_size = 100;
iterations = 50;
%父代选择的参数设置
parent_size = pop_size*0.6;
ud = 0.4; up = 1.6;
% 获取种群范围
[bounds,dimension] = get_variable_bounds(x);
%种群初始化
pop = init_pop(pop_size,dimension,bounds,x);
%种群排序
pop = sort_pop(pop,target,dimension);
% 初始化函数返回数据。
% my_nsga2_test.m 中需要保存的数据。
GD = zeros(1,iterations);
SP = zeros(1,iterations);
allpop = zeros(iterations,pop_size,dimension+4);%保存进化过程中种群的数据
warning off all
%迭代循环
for i = 1:iterations
%选择父代
parent_pop = select_parent(pop,ud,up,parent_size);
%进行遗传算法，杂交变异
child_pop = myga(parent_pop,dimension,bounds,i,iterations,x);
%子代和原始种群进行合并
pop = combined_pop(pop,child_pop,target,dimension);
%对合并种群进行非支配排序
pop = sort_pop(pop,target,dimension);
%选择新一代种群
pop = select_pop(pop,target,dimension,pop_size);
% 父代选择的参数变化
if ~mod(i,50)
if up < 2
ud = ud-0.05;
up = up+0.05;
end
end
%画出种群迭代的过程。只运行naga2_main的的时候，可以画出单个测试函数的变化
plot(pop(:,dimension+1),pop(:,dimension+2),'*')
grid on
title(['NSGA2测试第',num2str(x),'个函数第 ',num2str(i),' 代结果'])
pause(0.1)
%保存数据，计算每一代的GD和SP，也可以通过保存allpop后单独计算
allpop(i,:,:);
GD(1,i) = calculate_gd(pop,x);
SP(1,i) = calculate_sp(pop);
end
% 与理想前沿图画比较
plotPareto(x);
end
get_variable_bounds
function [bounds,dimension] = get_variable_bounds(x)
switch x
case 1
dimension = 30;
bounds = [ones(dimension,1)*0,ones(dimension,1)*1];
case 2
dimension = 30;
bounds = [ones(dimension,1)*0,ones(dimension,1)*1];
case 3
dimension = 30;
bounds = [ones(dimension,1)*0,ones(dimension,1)*1];
case 4
dimension = 10;
bounds = [zeros(1,1),ones(1,1);ones(9,1).*-5,ones(9,1).*5];
case 5
dimension = 10;
bounds = [ones(dimension,1)*0,ones(dimension,1)*1];
end
init_pop
function pop = init_pop(pop_size,dimension,bounds,x)
%佳点集生成初始种群
p = zeros(pop_size,dimension);
prime_number_min = dimension*2 +3;
% 找到(prime_number_min-3)/2>=dimension的最小素数prime_number_min
while 1
if isprime(prime_number_min)==1
break;
else
prime_number_min = prime_number_min + 1;
end
end
for i = 1:pop_size
for j = 1:dimension
r = mod(2*cos(2*pi*j/prime_number_min)*i,1);
p(i,j) = bounds(j,1)+r*(bounds(j,2)-bounds(j,1));
end
end
%计算种群的适应值
evaluate = calculate_pop(p,x);
pop = [p,evaluate];
sort_pop
function pop = sort_pop(pop_eva,target,dimension)
[N, ~] = size(pop_eva);
front = 1;
F(front).f = [];
individual = [];
%先确定等级为1的个体以及被支配的集合
for i = 1:N
individual(i).n = 0; %支配i的个体个数
individual(i).p = [];%被个体i支配的个体集合
for j = 1:N
less = 0; %判断i是否可以支配j
equal = 0; %判断i是否等于j，序号相同时相等
more = 0; %判断i是否被j支配
for k = 1:target %在每一个目标函数中判断支配关系
if pop_eva(i,dimension+k) < pop_eva(j,dimension+k)
less = less+1;
elseif pop_eva(i,dimension+k) == pop_eva(j,dimension+k)
equal = equal+1;
else
more = more + 1;
end
end
if less == 0 && equal ~= target
individual(i).n = individual(i).n + 1;
elseif more == 0 && equal ~= target
individual(i).p = [individual(i).p j];
end
end
if individual(i).n == 0
pop_eva(i,target+dimension+1) = 1;
F(front).f = [F(front).f i];
end
end
%对对所有种群所有个体进行等级划分
while ~isempty(F(front).f)
Q = [];
for i = 1:length(F(front).f) %等级为1的长度
if ~isempty(individual(F(front).f(i)).p) %等级为1的个体中查找其所支配的个体
for j = 1:length((individual(F(front).f(i)).p))%当前个体等级为1的个体所支配的个体数量
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
pop_eva(individual(F(front).f(i)).p(j),target + dimension + 1) = front + 1;
Q = [Q individual(F(front).f(i)).p(j)]; %记录下一等级的集合
end
end
end
end
front = front + 1;
F(front).f = Q;
end
%排序
[~, index_front] = sort(pop_eva(:,target + dimension +1));%根据等级对个体进行排序
sort_front = zeros(size(pop_eva));
for i = 1 : length(index_front)
sort_front(i,:) = pop_eva(index_front(i),:); %排序后的结果
end
current_index = 0; %当前下标。
%计算拥挤距离
for front = 1 : (length(F)-1)
distance = 0;
y =[];
previous_index = current_index + 1;
for i = 1 : length(F(front).f)
y(i,:) = sort_front(current_index + i,:);
end
current_index = current_index + i;
sorted_based_on_objective = [];
%函数值排序
for i = 1 : target
%函数值排序
[sorted_based_on_objective, index_of_objectives] = sort(y(:,dimension + 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), dimension + i);
f_min = sorted_based_on_objective(1, dimension + i);
y(index_of_objectives(length(index_of_objectives)),target + dimension + 1 + i)...
= Inf;
y(index_of_objectives(1),target + dimension + 1 + i) = Inf;
for j = 2 : length(index_of_objectives) - 1
next_obj = sorted_based_on_objective(j + 1,dimension + i);
previous_obj = sorted_based_on_objective(j - 1,dimension + i);
if (f_max - f_min == 0)
y(index_of_objectives(j),target + dimension + 1 + i) = Inf;
else
y(index_of_objectives(j),target + dimension + 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 : target
distance(:,1) = distance(:,1) + y(:,target + dimension + 1 + i);
end
y(:,target + dimension + 2) = distance;
y = y(:,1 : target + dimension + 2);
z(previous_index:current_index,:) = y;
end
pop = z();
select_pop
function pop = select_pop(pop,target,dimension,pop_size)
[popsize,~] = size(pop);
%根据等级对pop进行升序排序，对拥挤距离进行降序排序
sort_pop = sortrows(pop,[target+dimension+1,-(target+dimension+2)]);
s_pop = [];
no_index = [];
num = 0;
for i = 1:popsize-1
a = sort_pop(i,dimension+1:dimension+2);
b = sort_pop(i+1,dimension+1:dimension+2);
if norm(a-b)>1e-10 % 如果间隔大于1e-10，就选择当前个体
s_pop = [s_pop;sort_pop(i,:)];
num = num+1;
if num == pop_size
break;
end
else % 当前个体与下一个个体的间隔小于1e-10，就记录当前个体的索引
no_index = [no_index;i];
end
end
if size(s_pop,1)< pop_size
n = pop_size - size(s_pop,1);
for j = 1:n
s_pop = [s_pop;sort_pop(no_index(j),:)];
end
end
pop = s_pop;
select_parent
function parent_pop = select_parent(pop,ud,up,parent_size)
% 采用线性排名选择父代个体
[sizepop,td] = size(pop);
%对种群的等级和拥挤距离进行排序
pop = sortrows(pop,[-(td-1),td]);
parent_pop = [];
pi = zeros(1,sizepop);
for i = 1:sizepop
%计算每一个个体选择的概率
pi(1,i) = 1/sizepop*(ud+(up-ud)*((i-1)/(sizepop-1)));
end
pi = cumsum(pi);
r = sort(rand(1,parent_size));
fitin = 1; % 表示第几个个体
newin = 1; % 表示第几个随机数
while newin<=parent_size
if(r(newin)) <= pi(fitin)
parent_pop = [parent_pop;pop(fitin,:)];
newin = newin+1;
else
fitin = fitin+1;
end
end
myga
function child_pop = myga(parent_pop,dimension,bounds,t,T,x)
% t当前进进化次数，T最大迭代次数，x，问题编号
%GA算法
parent_pop = sortrows(parent_pop,[2+dimension+1,-(2+dimension+2)]);
parent_pop = parent_pop(:,1:dimension);
[popsize,~] = size(parent_pop);
%定义交叉变异的概率
crossover = 1;
mutation = 1;
%变异参数。b影响变异的范围。
b = 3;
child = [];
for i = 1:popsize
%判断个体是否进行遗传操作
c_r = rand(1);
m_r = rand(1);
%进行杂交操作
if c_r < crossover
%随机选择一个个体与该个体进行杂交
p1 = randperm(popsize,1);
if p1<=i
parent1 = parent_pop(p1,:);
parent2 = parent_pop(i,:);
else
parent1 = parent_pop(i,:);
parent2 = parent_pop(p1,:);
end
child1 = zeros(1,dimension);
child2 = zeros(1,dimension);
vga1 = parent1 - parent2;
vga2 = (parent1 + parent2)*0.5;
v1 = 0.729*vga1 + 1.5*rand*(parent1-vga2) + rand*(parent2 - vga2) ;
v2 = 0.729*vga1 + 1.5*rand*(parent2-vga2) + rand*(parent1 - vga2) ;
child1 = parent1 + v1;
child2 = parent2 + v2;
for j = 1:dimension
if child1(j) > bounds(j,2)
child1(j) = bounds(j,2);
elseif child1(j) < bounds(j,1)
child1(j) = bounds(j,1);
end
if child2(j) > bounds(j,2)
child2(j) = bounds(j,2);
elseif child2(j) < bounds(j,1)
child2(j) = bounds(j,1);
end
end
child = [child;child1;child2];
end
if m_r < mutation
% 非均匀变异
k = randperm(dimension,1);
rk = rand(1);
child3 = parent_pop(i,:);
if rk >= 0.5
y = bounds(k,2) - parent_pop(i,k);
xk = parent_pop(i,k) + y*(1-rand^((1-t/T)^b));
child3(1,k) = xk;
else
y = parent_pop(i,k) - bounds(k,1);
xk = parent_pop(i,k) - y*(1-rand^((1-t/T)^b));
child3(1,k) = xk;
end
child = [child;child3];
end
end
child_eva = calculate_pop(child,x);
child_pop = [child,child_eva];
combined_pop
function pop = combined_pop(pop,child_pop,target,dimension)
%合并父代和子代个体
pop1 = pop(:,1:target+dimension);
clear pop
pop = [pop1;child_pop];
calculate_gd
function GD = calculate_gd(pop,x)
switch x
case 1
y = importdata('前沿数据/ZDT1.txt');
case 2
y = importdata('前沿数据/ZDT2.txt');
case 3
y = importdata('前沿数据/ZDT3.txt');
case 4
zdt4 = importdata('前沿数据/ZDT4.txt');
y = sortrows(zdt4,[1,2]);
case 5
y = importdata('前沿数据/ZDT6.txt');
end
%pop测试结果，y真实值
GD = 0;
[n,d] = size(pop);
pop = pop(:,d-3:d-2);
for i = 1:n
dis = pdist2(pop(i,:),y,'euclidean');
gd = (min(dis))^2;
% gd = min(dis);
GD = GD + gd;
end
GD = sqrt(GD/n);
% GD = GD/n
end
calculate_sp
function SP = calculate_sp(pop)
[x,y] = size(pop);
pop = pop(:,y-3:y-2);
mindis = zeros(x,1);
for i = 1:x
di = pop(i,:);
dis = pdist2(di,pop,'euclidean');
dis = sort(dis);
mindis(i) = dis(2);
end
meandis = mean(mindis);
Sp = 0;
for j = 1:x
sp = (meandis-mindis(j))^2;
Sp = Sp + sp;
end
SP = sqrt(Sp/x)/meandis;
calculate_pop
function evaluate = calculate_pop(pop,x)
%测试函数
%global x
[~,dim] = size(pop);
switch x
case 1 %ZDT1
fx1 = pop(:,1);
gx = 1+sum(pop(:,2:end),2).*(9/(dim-1));
hx = 1-sqrt(fx1./gx);
fx2 = gx.*hx;
evaluate = [fx1,fx2];
case 2 %ZDT2
fx1 = pop(:,1);
gx = 1+sum(pop(:,2:end),2).*(9/(dim-1));
hx = 1-(fx1./gx).^2;
fx2 = gx.*hx;
evaluate = [fx1,fx2];
case 3 %ZDT3
fx1 = pop(:,1);
gx = 1+sum(pop(:,2:end),2).*(9/(dim-1));
hx = 1-sqrt(fx1./gx)-(fx1./gx).*sin(10*pi.*fx1);
fx2 = gx.*hx;
evaluate = [fx1,fx2];
case 4 %ZDT4
fx1 = pop(:,1);
gx = 91+sum((pop(:,2:dim).^2-10.*cos(4*pi.*pop(:,2:dim))),2);
hx = 1-sqrt(fx1./gx);
fx2 = gx.*hx;
evaluate = [fx1,fx2];
case 5
x1 = pop(:,1);
fx1 = 1-exp(-4.*x1).*(sin(6*pi.*x1)).^6;
s = sum(pop(:,2:end),2);
gx = 1+9/(dim-1).*s;
hx = 1-(fx1./gx).^2;
fx2 = gx.*hx;
evaluate = [fx1,fx2];
end
plotPareto
function plotPareto(x)
switch x
case 1
zdt1 = importdata('前沿数据/ZDT1.txt');
hold on
plot(zdt1(:,1),zdt1(:,2),'-')
legend('改进NSGA2测试前沿','理想前沿')
case 2
zdt2 = importdata('前沿数据/ZDT2.txt');
hold on
plot(zdt2(:,1),zdt2(:,2),'-')
legend('测试前沿','已知前沿')
case 3
zdt3 = importdata('前沿数据/ZDT3.txt');
hold on
plot(zdt3(:,1),zdt3(:,2),'*')
legend('测试前沿','已知前沿')
case 4
zdt4 = importdata('前沿数据/ZDT4.txt');
zdt4 = sortrows(zdt4,[1,2]);
hold on
plot(zdt4(:,1),zdt4(:,2),'-')
legend('测试前沿','已知前沿')
case 5
zdt6 = importdata('前沿数据/ZDT6.txt');
hold on
plot(zdt6(:,1),zdt6(:,2),'-')
legend('测试前沿','已知前沿')
otherwise
fprintf('错误')
end
end
运行结果
my_nsga2_main测试
测试zdt1，zdt2，zdt3，zdt4和zdt6的结果
参数设置为：
- pop_size = 100;
- iterations = 50;
-ZDT4函数测试比较特殊，需要的计算比较大，要得到靠进前沿的结果， 参数需要改为：pop_size =200；terations = 300;
ZDT1  | ZDT2 |
|---|---|
ZDT3 | ZDT4  |
ZDT6  | |
my_nsga2_test运行结果
注意事项：
统一参数设置：由于zdt4函数的特殊，建议统一改为
- pop_size = 100;
- iterations = 50;
my_nsga2_test 中的参数设置，与my_nsga2_main的只需要设置一处即可。
此脚本主要是为了保存测试过程产生的数据。不需要查看种群迭代的情况，将画图相关的可以注释掉，同时建议测试次数test改为多次，如30。
个人分享过程就只设置为1了。
 |  |
|---|
1187
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
