多目标优化NSGA-II的实现（MATLAB完整代码）_matlab nsga2

由于历史原因，没有整理好完整的代码，所以在【多目标优化NSGA-II的实现和测试（MATLAB实现）】中只放了部分代码。

现在已经整理好了代码，此部分的代码测试内容为：ZDT1、ZDT2、ZDT3、ZDT4、ZDT6。

---

目录

主要内容

代码模块

其他内容

运行注意事项

 代码

nsga2_test

nsga2_main

get_variable_bounds

init_pop

sort_pop

select_parent

myga

combined_pop

select_pop

calculate_gd

calculate_sp

calculate_pop

plotPareto

运行结果

---

---

主要内容

代码模块

• nsga2_test：测试函数，用于保存测试数据
• 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-互联网文档类资源-CSDN文库https://download.csdn.net/download/weixin_44034444/73514580

运行注意事项

1. 在nsga2_test中设置pop_size，iterations。以及测试次数test。nsga2_test可以保存每一个测试函数，每一测试中每一代的种群数据以及GD和SP的数据。主要是为了方便用获取得到的数据进行分析。
2. 如果只是需要查看nsga2的效果，运行nsga2_main函数，注意pop_size，iterations的设置。

 代码

nsga2_test

clc
clear
 
% 定义全局变量,
global pop_size; 
global iterations;
pop_size = 100;%种群大小
iterations = 500;%迭代次数
 
test = 1; %测试次数
 
for x = 1:5%选择需要计算的函数
    switch x
        case 1
            [~,dim] = get_variable_bounds(x);
            %dim = 10;
            %设置保存参数
            testGD = zeros(test,iterations);
            testSP = zeros(test,iterations);
            testPop = zeros(test,iterations,pop_size,dim+4);
            NSGA2_zdt1 = [];
            disp('正在测试zdt1')
            for i = 1:test
                disp(['第' num2str(i) '次测试']);
                [pop,GD,SP] = nsga2_main(x);
                testGD(i,:) = GD;
                testSP(i,:) = SP;
                testPop(i,:,:,:) = pop;      
            end
            NSGA2_zdt1.testGD = testGD;
            NSGA2_zdt1.testSP = testSP;
            NSGA2_zdt1.testPop = testPop;
            save('NSGA2_zdt1.mat','NSGA2_zdt1')
        case 2
            [~,dim] = get_variable_bounds(x);
            %dim = 10;
            %设置保存参数
            testGD = zeros(test,iterations);
            testSP = zeros(test,iterations);
            testPop = zeros(test,iterations,pop_size,dim+4);
            NSGA2_zdt2 = [];
            disp('正在测试zdt2')
            for i = 1:test
                disp(['第' num2str(i) '次测试']);
                [pop,GD,SP] = nsga2_main(x);
                testGD(i,:) = GD;
                testSP(i,:) = SP;
                testPop(i,:,:,:) = pop;      
            end
            NSGA2_zdt2.testGD = testGD;
            NSGA2_zdt2.testSP = testSP;
            NSGA2_zdt2.testPop = testPop;
            save('NSGA2_zdt2.mat','NSGA2_zdt2')
       case 3 
            [~,dim] = get_variable_bounds(x);
            %dim = 10;
            %设置保存参数
            testGD = zeros(test,iterations);
            testSP = zeros(test,iterations);
            testPop = zeros(test,iterations,pop_size,dim+4);
            NSGA2_zdt3 = [];
            disp('正在测试zdt3')
            for i = 1:test
                disp(['第' num2str(i) '次测试']);
                [pop,GD,SP] = nsga2_main(x);
                testGD(i,:) = GD;
                testSP(i,:) = SP;
                testPop(i,:,:,:) = pop;      
            end
            NSGA2_zdt3.testGD = testGD;
            NSGA2_zdt3.testSP = testSP;
            NSGA2_zdt3.testPop = testPop;
            save('NSGA2_zdt3.mat','NSGA2_zdt3')
        case 4
            [~,dim] = get_variable_bounds(x);
            %dim = 10;
            %设置保存参数
            testGD = zeros(test,iterations);
            testSP = zeros(test,iterations);
            testPop = zeros(test,iterations,pop_size,dim+4);
            NSGA2_zdt4 = [];
            disp('正在测试zdt4')
            for i = 1:test
                disp(['第' num2str(i) '次测试']);
                [pop,GD,SP] = nsga2_main(x);
                testGD(i,:) = GD;
                testSP(i,:) = SP;
                testPop(i,:,:,:) = pop;      
            end
            NSGA2_zdt4.testGD = testGD;
            NSGA2_zdt4.testSP = testSP;
            NSGA2_zdt4.testPop = testPop;
            save('NSGA2_zdt4.mat','NSGA2_zdt4')
        case 5
            [~,dim] = get_variable_bounds(x);
            %dim = 10;
            %设置保存参数
            testGD = zeros(test,iterations);
            testSP = zeros(test,iterations);
            testPop = zeros(test,iterations,pop_size,dim+4);
            NSGA2_zdt6 = [];
            disp('正在测试zdt6')
            for i = 1:test
                disp(['第' num2str(i) '次测试']);
                [pop,GD,SP] = nsga2_main(x);
                testGD(i,:) = GD;
                testSP(i,:) = SP;
                testPop(i,:,:,:) = pop;      
            end
            NSGA2_zdt6.testGD = testGD;
            NSGA2_zdt6.testSP = testSP;
            NSGA2_zdt6.testPop = testPop;
            save('NSGA2_zdt6.mat','NSGA2_zdt6')
    end
    
end

nsga2_main

function [allpop,GD,SP] = nsga2_main(x)
% 测试主函数 x,问题编号
% 输出种群，GD和SP
 
% 参数设置 
global pop_size
global iterations;%迭代次数
target = 2;
 
% 获取种群范围
[bounds,dimension] = get_variable_bounds(x);
%种群初始化
pop = init_pop(pop_size,dimension,bounds,x);
%种群排序
pop = sort_pop(pop,target,dimension);
 
%锦标赛参数设置
parent_size = pop_size/2;
select_size = 2;
 
% 初始化函数返回数据。
% nsga2_test.m 中需要保存的数据。 如果不跑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,parent_size,select_size);
    %进行遗传算法，杂交变异
    child_pop = myga(parent_pop,dimension,bounds,x);
    %子代和原始种群进行合并
    pop = combined_pop(pop,child_pop,target,dimension);
    %对合并种群进行非支配排序
    pop = sort_pop(pop,target,dimension);
    %选择新一代种群
    pop = select_pop(pop,target,dimension,pop_size);
    
%    %画出种群迭代的过程。只运行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,:,:) = pop;
    GD(1,i) = calculate_gd(pop,x);
    SP(1,i) = calculate_sp(pop);
end
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 = rand(pop_size,dimension);%生成popsize*dimension的0-1矩阵
%生成定义域范围内种群
for i = 1:dimension
    p(:,i) = bounds(i,1)+p(:,i)*(bounds(i,2)-bounds(i,1));
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_parent

function parent_pop = select_parent(pop,parent_size,compare_size)
%父代个体的选择
[pop_size,distance] = size(pop);
rank = distance-1; %记录等级所在的列
select_pop = zeros(compare_size,distance);
 
for i = 1:parent_size
    %生成参与锦标赛的个体序列
    parent_list = randperm(pop_size,compare_size);
    %参与锦标赛的个体集合
    for j = 1:compare_size
        select_pop(j,:) = pop(parent_list(j),:);
    end
    [min_rank,min_rank_index] = min(select_pop(:,rank));
    if length(min_rank)==1
        parent_pop(i,:) = select_pop(min_rank_index,:);
    else
        %最小等级相同的个体集合
        for k = 1:length(min_rank)
            select_pop1(k,:) = select_pop(min_rank_index(k),:);
        end
        [~,max_distance_index] = max(select_pop1(:,distance));
        parent_pop(i,:) = select_pop1(max_distance_index(1),:);
    end
end

myga

function child_pop = myga(parent_pop,dimension,bounds,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;
nc=20;
child = [];
 
 
for i = 1:popsize
    c_r = rand(1);
    m_r = rand(1);
    %交叉变换
    if c_r < crossover
        %随机选择一个个体与该个体进行杂交
        p1 = randperm(popsize,1); 
        parent1 = parent_pop(p1,:);
        parent2 = parent_pop(i,:);
           
        % 多项式杂交
        child1 = zeros(1,dimension);
        child2 = zeros(1,dimension);
        for j = 1:dimension
            r = rand(1);
            if r <= 0.5
                a = (2*r)^(1/(nc+1));
            else
                a= (2*(1-r))^(-(1/(nc+1)));
            end
            
            child1(j) = ((1+a)*parent1(j) + (1-a)*parent2(j))/2;
            child2(j) = ((1-a)*parent1(j) + (1+a)*parent2(j))/2;
            
            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
        child3=parent_pop(i,:);
        for k = 1:dimension
            r = rand();
            if r<0.5
                m = (2*r)^(1/21)-1;
            else
                m = 1 - (2*(1 - r))^(1/(21));
            end
            child3(1,k) = child3(1,k)+m;
            if child3(1,k)>bounds(k,2)
                child3(1,k) = bounds(k,2);
            end
            if child3(1,k)<bounds(k,1)
                child3(1,k) = bounds(k,1);
            end
        end  
        child = [child;child3];
    end
end
child_pop = child(:,1:dimension);
child_eva = calculate_pop(child_pop,x);
child_pop = [child_pop,child_eva];

combined_pop

function pop = combined_pop(pop,child_pop,target,dimension)
%合并父代和子代个体
pop1 = pop(:,1:target+dimension);
clear pop
pop = [pop1;child_pop];

select_pop

function pop = select_pop(pop,target,dimension,pop_size)
[popsize,~] = size(pop);
sort_pop = sortrows(pop,[target+dimension+1,-(target+dimension+2)]);
s_pop = [];
no_index = [];
num = 0;
if popsize > pop_size
    %根据等级对pop进行升序排序，对拥挤距离进行降序排序   
    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
            s_pop = [s_pop;sort_pop(i,:)];
            num = num+1;
            if num == pop_size
                break;
            end
        else
            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;
end

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)
%测试函数
[~,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];
    case 6
        n = -sum((pop-1/sqrt(dim)).^2,2);
        m = -sum((pop+1/sqrt(dim)).^2,2);
        fx1 = 1-exp(n);
        fx2 = 1-exp(m);
        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

运行结果

运行过程

保存的数据