多目标zdt测试问题—moead、nsga2、spea2算法

big_f.m

function [r] = big_f(i,j)
    global rank;
    global distance;
    if rank(i) > rank(j) || (rank(i) == rank(j) && distance(i) < distance(j) )
        r = true;
    else
        r = false;
    end
end


bip.m

function [r] = bip(f,lamda,z,M,theta)
    d1 = norm((f-z)'*lamda)/norm(lamda);
    d2 = norm(f-(z+d1*lamda));
    r = d1 + theta * d2;
end


crowding_distance_assignment.m

function [] = crowding_distance_assignment(I,f)
    l = size(I,2);
    if l < 1
        return;
    end
    global distance;
    for m = 1:2
        [b,idx] = sort(f(m,:));
        distance(1,idx(1)) = Inf;
        distance(1,idx(l)) = Inf;
        for i = 2:(l-1)
          distance(1,idx(i)) = distance(1,idx(i)) + (b(1,i+1) - b(1,i-1)) / (b(1,l) - b(1,1));
        end
    end
end


distance_relation.m

function [r] = distance_relation(Arg1,Arg2)
    for i = 1:size(Arg1,2)
        if Arg1(i) > Arg2(i)
            r = -1;
            return;
        end
        if Arg1(i) < Arg2(i)
            r = 1;
            return;
        end
    end
    r = 0;
end


dominate.m

function [r] = dominate(x,y,M)
    lte = 0;
    lt = 0;
    gte = 0;
    gt = 0;
    for i = 1:M
        if x(i)<=y(i) 
            lte=lte+1; 
        end
        if x(i)<y(i) 
            lt=lt+1; 
        end
        if x(i)>=y(i) 
            gte=gte+1; 
        end
        if x(i)>y(i) 
            gt=gt+1; 
        end
    end
    if lte==M && lt>0
        r=1;
    elseif gte==M && gt>0
        r=-1;
    else
        r=0;
    end
end


dominate_relation.m

function [r] = dominate_relation(x,y,f)
    num_gt = 0;
    num_gte = 0;
    num_lt = 0;
    num_lte = 0;
    for i = 1:2
        if f(i,x) < f(i,y) num_lt=num_lt+1; end
        if f(i,x) <= f(i,y) num_lte=num_lte+1; end
        if f(i,x) > f(i,y) num_gt=num_gt+1; end
        if f(i,x) >= f(i,y) num_gte=num_gte+1; end
    end
    if num_lte == 2 && num_lt > 0
        r = -1;%支配
        return;
    end
    if num_gte == 2 && num_gt > 0
        r = 1;%被支配
        return;
    end
    r = 0;
end


fast_non_dominated_sort.m

function [] = fast_non_dominated_sort(P,f)
    global rank;
    rank = [];
    global F;
    for p = 1:size(P,1)
        s(p).set = [];
        n(p) = 0;
        for q = 1:size(P,1)
            r = dominate_relation(p,q,f);
            if r == -1
               s(p).set = [s(p).set q];
            elseif r == 1
               n(p) = n(p) + 1;
            end
        end
        if n(p) == 0
            rank(p) = 1;
            F(1).set = [F(1).set p];
        end
    end
    i = 1;
    while size(F(i).set,2) > 0
        Q = [];
        for ii = 1:size(F(i).set,2)
            p = F(i).set(1,ii);
            for j = 1:size(s(p).set,2)
                q = s(p).set(1,j);
                n(q) = n(q) - 1;
                if n(q) == 0
                    rank(q) = i + 1;
                    Q = [Q q];
                end
            end
        end
        i = i + 1;
        F(i).set = Q; 
    end
end


gte.m

function [r] = gte(f,lamda,z,M)
    r = max(lamda.*abs(f-z));
end


make_new_pop.m

function [q] = make_new_pop(p,num_geti,num_x,type,u_limit,d_limit)
    size_p = size(p,1);
    mu = 20;
    idx_now = 1;
    for half_idx = 1:fix(num_geti/2)
       % if rand(1) < 0.9
        father_idx = select_one_bytour(size_p,type);
        mother_idx = select_one_bytour(size_p,type);
        while father_idx == mother_idx
             mother_idx = select_one_bytour(size_p,type);
        end
        son_1 = idx_now;
        son_2 = idx_now+1;
        for j = 1:num_x
            u = rand(1);
            if u < 0.5
                yita = (2*u)^(1/(mu+1));
            else
                yita = (1/(2-2*u))^(1/(mu+1));
            end
            q(son_1,j) = 0.5*((1+yita)*p(father_idx,j)+(1-yita)*p(mother_idx,j));
            q(son_2,j) = 0.5*((1-yita)*p(father_idx,j)+(1+yita)*p(mother_idx,j)); 
            if rand(1) < 0.2
                if u < 0.5
                    delta = (2*u)^(1/(mu+1)) - 1;
                else
                    delta = 1-(2-2*u)^(1/(mu+1));
                end
                delta = delta*(u_limit(1,j)-d_limit(1,j));
                q(son_1,j) = q(son_1,j) + delta;
            end
            if rand(1) < 0.2
                if u < 0.5
                    delta = (2*u)^(1/(mu+1)) - 1;
                else
                    delta = 1-(2-2*u)^(1/(mu+1));
                end
                delta = delta*(u_limit(1,j)-d_limit(1,j));
                q(son_2,j) = q(son_2,j) + delta;
            end
            if q(son_1,j) < d_limit(1,j)
                q(son_1,j) = d_limit(1,j);
            end
            if q(son_2,j) < d_limit(1,j)
                q(son_2,j) = d_limit(1,j);
            end
            if q(son_1,j) > u_limit(1,j)
                q(son_1,j) = u_limit(1,j);
            end
            if q(son_2,j) > u_limit(1,j)
                q(son_2,j) = u_limit(1,j);
            end     
        end
        idx_now = idx_now + 2;
       % else
%             father_idx = select_one_bytour(size_p,type);
%             son_1 = idx_now;
%             for j = 1:num_x
%                 u = rand(1);
%                 if u < 0.5
%                     delta = (2*u)^(1/(mu+1)) - 1;
%                 else
%                     delta = 1-(2-2*u)^(1/(mu+1));
%                 end
%                 delta = delta*(u_limit(1,j)-d_limit(1,j));
%                 q(son_1,j) = p(father_idx,j) + delta;
%                 if q(son_1,j) < d_limit(1,j)
%                     q(son_1,j) = d_limit(1,j);
%                 end
%                 if q(son_1,j) > u_limit(1,j)
%                     q(son_1,j) = u_limit(1,j);
%                 end
%             end
%             idx_now = idx_now + 1;
%         end
    end
end


MOEAD主程序moead.m

%cd('../');
%addpath(genpath(cd));
clear;
num_shidai = 500;
N = 100;
M = 2;
lamda = zeros(N,M);
for i = 1:N
    lamda(i,1) = i/N;
    lamda(i,2) = (N-i)/N;
end
T = fix(N/5);
num_x = 30;
mu = 20;
theta = 5;%for 

EP = [];
for i = 1:N
    for j = 1:N
        B(i,j) = norm(lamda(i,:)-lamda(j,:));
        B(j,i) = B(i,j);
    end
end
[~,B] = sort(B,2);%按列排序

now_number = 0;
for zdttype = 2:2
    for method = 0:2 %0:Weight Sum Approach  1:PBI 2:chebycheff
        if 5 == zdttype
            continue;
        end
        x=[];
        EP=[];
        if zdttype >= 4
            num_x = 10;
        else
            num_x = 30;
        end
        if zdttype~=4
            x = rand(N,num_x);
        else
            x = zeros(N,num_x);
            x(:,1) = rand(N,1);
            x(:,2:num_x) = rand(N,num_x-1)*10-5;
        end
        for i = 1:N
            [f(1,i),f(2,i)] = zdt(x(i,:),zdttype);
        end
        for i = 1:M
           z(i) = min(f(i,:));
           %z(i) = max(f(i,:));
        end
        for i = 1:N
            flag = 1;
            for j = 1:N
                if 1 == dominate(f(:,j),f(:,i),M)
                    flag = 0;
                    break;
                end
            end
            if flag == 1
                EP = [EP; x(i,:) f(:,i)'];
            end
        end
        for shidailoop = 1:num_shidai
            for i = 1:N
                k = B(i,randi(T));
                l = B(i,randi(T));
                while k == l
                    l = B(i,randi(T));
                end
                y = zeros(1,num_x);
                for j = 1:num_x  
                    u = rand(1);
                    if u < 0.5
                        yita = (2*u)^(1/(mu+1));
                    else
                        yita = (1/(2-2*u))^(1/(mu+1));
                    end
                    if rand(1) < 0.5
                        y(1,j) = 0.5*((1+yita)*x(k,j)+(1-yita)*x(l,j));
                    else
                        y(1,j) = 0.5*((1-yita)*x(k,j)+(1+yita)*x(l,j));
                    end
                    d_limit = 0;
                    u_limit = 1;
                    if zdttype == 4 && j>1
                        d_limit = -5;
                        u_limit = 5;
                    end
                    if rand(1) < 0.2
                        u = rand(1);
                        if u < 0.5
                            delta = (2*u)^(1/(mu+1)) - 1;
                        else
                            delta = 1-(2-2*u)^(1/(mu+1));
                        end
                       y(1,j) = y(1,j) + delta*(u_limit-d_limit);
                    end
                    if y(1,j) < d_limit
                        y(1,j) = d_limit;
                    end
                    if y(1,j) > u_limit
                        y(1,j) = u_limit;
                    end
                end
                [fj(1),fj(2)] = zdt(y,zdttype);
                for j = 1:M
                    if z(j) > fj(j)
                        z(j) = fj(j);
                    end
                end
                for j = 1:T
                    p = B(i,j);
                    if 2 == method
                       value_fj = gte(fj,lamda(p,:),z,M);
                       value_p = gte(f(:,p)',lamda(p,:),z,M);
                    elseif 1 == method
                       value_fj = bip(fj,lamda(p,:),z,M,theta);
                       value_p = bip(f(:,p)',lamda(p,:),z,M,theta);
                    elseif 0 == method
                       value_fj = sum(fj.*lamda(p,:));
                       value_p = sum(f(:,p)'.*lamda(p,:));
                    end
                    if value_fj < value_p
                    %if value_fj > value_p
                        x(p,:) = y;
                        for loop = 1:M
                            f(loop,p) = fj(loop);
                        end
                    end
                end
                flag=1;
                for j = 1:size(EP,1)
                    if j>size(EP,1)
                        break;
                    end
                    r = dominate(fj,EP(j,num_x+1:end),M);
                    if 1 == r
                   % if -1 == r
                        EP(j,:) = [];
                    elseif -1 == r
                   % elseif 1 == r
                        flag=0;
                    end
                end
                if flag==1
                    EP = [EP;y fj];
                end
            end
            figure(method*6+zdttype);
            cla;
            plot(EP(:,num_x+1),EP(:,num_x+2),'*');
            xlabel('f_1'); ylabel('f_2');
            hold on;
            optarg = zeros(100,num_x);
            optarg(:,1) =  linspace(0,1);
            xx = zeros(1,100);
            yy = zeros(1,100);
            for i = 1:100
                [xx(1,i) yy(1,i)]=zdt(optarg(i,:),zdttype);
            end
            plot(xx,yy);
            if method == 2
                str = ['ZDT', num2str(zdttype),'-Chebycheff'];
            elseif method == 1
                str = ['ZDT', num2str(zdttype),'-PBI'];
            elseif method == 0
                str = ['ZDT', num2str(zdttype),'-Weight Sum Approach'];
            end
            title(str);
            pause(0.01);
            if ~mod(shidailoop,50)
                clc;
                fprintf('%d generations completed\n',shidailoop);
            end
        end
        shidailoop = shidailoop;
    end
end

NSGA2主程序nsga2.m

num_geti = 100;
num_x = 30;
num_shidai = 500;

global p_t;
global q_t;
global r_t;
global F;
global rank;
global distance;

for zdttype = 1:6
    if 5 == zdttype
        continue;
    end
    if zdttype >= 4
        num_x = 10;
    end
    if zdttype~=4
        p_t = rand(num_geti,num_x);
    else
        p_t = zeros(num_geti,num_x);
        p_t(:,1) = rand(num_geti,1);
        p_t(:,2:num_x) = rand(num_geti,num_x-1)*10-5;
    end
    q_t = [];
    u_limit = ones(1,num_x);
    d_limit = zeros(1,num_x);
    if (4 == zdttype)
        for i = 2:num_x
            u_limit(i) = 5;
            d_limit(i) = -5;
        end
    end
    for shidai_loop = 1:num_shidai    
        r_t = [p_t;q_t];
        sum_size = size(r_t,1);
        f = zeros(2,sum_size);
        for i = 1:sum_size
            [f(1,i),f(2,i)] = zdt(r_t(i,:),zdttype);
        end
        for i = 1:sum_size
            F(i).set = [];
        end
        fast_non_dominated_sort(r_t,f);
        p_t1_idx = []; 
        i = 1;
        while size(p_t1_idx,2) + size(F(i).set,2) <= num_geti && size(F(i).set,2) > 0
            p_t1_idx = [p_t1_idx F(i).set];
            i = i+1;
        end
        for idx = 1:sum_size
            distance(1,idx) = 0;
        end
       crowding_distance_assignment(F(i).set,f);
       sort_f(i,1,size(F(i).set,2));
       p_t1_idx = [p_t1_idx F(i).set(1:num_geti-size(p_t1_idx,2))];
       p_t1 = r_t(p_t1_idx,:);
       q_t1 = make_new_pop(p_t1,num_geti,num_x,1,u_limit,d_limit);
       p_t = p_t1;
       q_t = q_t1;
       if ~mod(shidai_loop,50)
            clc;
            fprintf('%d generations completed\n',shidai_loop);
       end
       figure(zdttype);
       cla;
        size_ff = size(F(1).set,2);
        ff = zeros(2,size_ff);
        for i = 1:size_ff
            [ff(1,i),ff(2,i)] = zdt(r_t(F(1).set(i),:),zdttype);
        end
        plot(ff(1,:),ff(2,:),'*');
        xlabel('f_1'); ylabel('f_2');
        hold on;
        optarg = zeros(100,num_x);
        optarg(:,1) =  linspace(0,1);
        x = zeros(1,100);
        y = zeros(1,100);
        for i = 1:100
            [x(1,i) y(1,i)]=zdt(optarg(i,:),zdttype);
        end
        plot(x(1,:),y(1,:));
        str = ['ZDT', num2str(zdttype)];
        title(str);
        pause(0.01);
    end
end

select_one_bytour.m

function [idx] = select_one_bytour(size_p,type)
     x = unidrnd(size_p);
     y = unidrnd(size_p);
     while x == y
        y = unidrnd(size_p);
     end
     if type ~= 1
         idx = min(x,y);
         return;
     end
     global rank;
     if rank(x) < rank(y)
         idx = x;
         return;
     end
     if rank(x) > rank(y)
         idx = y;
         return;
     end
     global distance;
     if distance(x) > distance(y)
         idx = x;
         return;
     end
     idx = y;
end


small_f.m

function [r] = small_f(i,j)
    global rank;
    global distance;
    if rank(i) < rank(j) || (rank(i) == rank(j) && distance(i) > distance(j) )
        r = true;
    else
        r = false;
    end
end


sort_distance.m

function [] = sort_distance(l,r)
    global distance_change;
    i=l;
    j=r;
    mid = distance_change(fix((l+r)/2),:);
    while i<=j
        while -1 == distance_relation(distance_change(i,:),mid) 
            i = i+1; 
        end
        while 1 == distance_relation(distance_change(j,:),mid)
            j = j-1;
        end
        if i<=j
            t = distance_change(i,:);
            distance_change(i,:) = distance_change(j,:);
            distance_change(j,:) = t;
            i = i+1;
            j = j-1;
        end
    end
    if (l<j)
        sort_distance(l,j);
    end
    if (i<r)
        sort_distance(i,r);
    end
end


sort_f.m

function [s_out] = sort_f(idx,l,r)
    %s_out = s_in;
    global F;
    i = l;
    j = r;
    if l > j
        return;
    end
    mi = F(idx).set(1,fix((i+j)/2));
    while (i<=j)
        while small_f(F(idx).set(i),mi) 
            i = i + 1;
        end
        while big_f(F(idx).set(j),mi) 
            j = j - 1;
        end
        if i <= j
            t = F(idx).set(i);
            F(idx).set(i) = F(idx).set(j);
            F(idx).set(j) = t;
            i = i + 1;
            j = j - 1;
        end
    end
    if (l < j)
        sort_f(idx,l,j);
    end
    if (i < r)
        sort_f(idx,i,r);
    end
end


SPEA2主程序spea2.m

%input
population_size = 100;
archive_size = 100;
max_generations = 500;
num_x = 30;
M = 2;
k = round(sqrt(population_size + archive_size));
%output
A=[];
%initalization
for zdttype = 1:6
    if 5 == zdttype
        continue;
    end
    if zdttype >= 4
        num_x = 10;
    end
    if zdttype~=4
        pt=rand(population_size,num_x);
    else
        pt = zeros(population_size,num_x);
        pt(:,1) = rand(population_size,1);
        pt(:,2:num_x) = rand(population_size,num_x-1)*10-5;
    end
    pt_ba=[];
    u_limit = ones(1,num_x);
    d_limit = zeros(1,num_x);
    if (4 == zdttype)
        for i = 2:num_x
            u_limit(i) = 5;
            d_limit(i) = -5;
        end
    end
    for now_generation = 0:max_generations
        rt = [pt;pt_ba];%pt+pt_ba
        size_rt = size(rt,1);
        for i=1:size_rt
            [f1,f2] = zdt(rt(i,:),zdttype);
            rt(i,num_x+1:num_x+6) = [f1 f2 0 0 0 0];
        end
        %计算S
        for i=1:size_rt-1
            for j=i+1:size_rt
                r = dominate(rt(i,num_x+1:num_x+2),rt(j,num_x+1:num_x+2),M);
                if 1 == r           %i支配j
                    rt(i,num_x+3) = rt(i,num_x+3) + 1;
                elseif -1 == r           %j支配i
                    rt(j,num_x+3) = rt(j,num_x+3) + 1;
                end
            end
        end
        %计算R
        for i=1:size_rt-1
            for j=i+1:size_rt
                r = dominate(rt(i,num_x+1:num_x+2),rt(j,num_x+1:num_x+2),M);
                if -1 == r                 %j支配i
                    rt(i,num_x+4) = rt(i,num_x+4) + rt(j,num_x+3);
                elseif 1 == r                %i支配j
                    rt(j,num_x+4) = rt(j,num_x+4) + rt(i,num_x+3);
                end
            end
        end
        %计算任意两点间距离
        distance = zeros(size_rt,size_rt);
        for i=1:size_rt-1
            for j=i+1:size_rt
                distance(i,j) = norm(rt(i,num_x+1:num_x+2)-rt(j,num_x+1:num_x+2));
                distance(j,i) = distance(i,j);
            end
        end
        %计算D,F
        for i=1:size_rt
            [distance(i,:),~] = sort(distance(i,:));
            rt(i,num_x+5) = 1/(distance(i,5)+2);
            rt(i,num_x+6) = rt(i,num_x+4) + rt(i,num_x+5);
        end
        pareto_index = find(rt(:,num_x+6)<1);
        p_t1ba = rt(pareto_index,:);
        size_p_t1 = size(p_t1ba,1);
        if size_p_t1 < archive_size
            [rt,pareto_index] = sortrows(rt,num_x+6);
            p_t1ba = [p_t1ba;rt(size_p_t1+1:archive_size,:)];
        elseif size_p_t1 > archive_size
            distance = distance(pareto_index,:);
            for i = 1:size_p_t1 %  size_p_t1=size(distance,1)
                distance(i,size_rt+1) = i;
            end
            global distance_change;
            distance_change = distance;
            sort_distance(1,size_p_t1);
            p_t1ba = p_t1ba(fix(distance_change(1:archive_size,size_rt+1)),:);
        end
        figure(zdttype);
        cla;
        plot(p_t1ba(:,num_x+1),p_t1ba(:,num_x+2),'*');
        xlabel('f_1'); ylabel('f_2');
        if now_generation < max_generations
            p_t1 = make_new_pop(p_t1ba,archive_size,num_x,2,u_limit,d_limit);
            pt = p_t1(:,1:num_x);
            pt_ba = p_t1ba(:,1:num_x);
            if ~mod(now_generation,50)
                clc;
                fprintf('%d generations completed\n',now_generation);
            end
        end
        hold on;
        optarg = zeros(100,num_x);
        optarg(:,1) =  linspace(0,1);
        x = zeros(1,100);
        y = zeros(1,100);
        for i = 1:100
            [x(1,i) y(1,i)]=zdt(optarg(i,:),zdttype);
        end
        plot(x(1,:),y(1,:));
        str = ['ZDT', num2str(zdttype)];
        title(str);
        pause(0.01);
    end
end

zdt.m

function [f1,f2] = zdt(arg,zdttype)
    switch zdttype
        case 6
            f1 = 1-exp(-4*arg(1,1))*(sin(6*pi*arg(1,1)))^6;
        otherwise
            f1 = arg(1,1);
    end
    g = 0; 
    if zdttype>=4
        n = 10;
    else
        n = 30;
    end
    switch zdttype
        case 4
            for i = 2:n
                g = g+arg(1,i)^2-10*cos(4*pi*arg(1,i));
            end
            g = g + 1 + 10*(n-1);
        case 6
            for i = 2:n
               g = g+arg(1,i);
            end
            g = 1 + 9*(g/(n-1))^0.25;
        otherwise
            for i = 2:n
               g = g+arg(1,i);
            end
            g = 1 + g*9/(n-1);
    end
    switch zdttype
        case {1,4}
            f2 = g*(1-(arg(1,1)/g)^0.5);
        case 2
            f2 = g*(1-(arg(1,1)/g)^2);
        case 3
            f2 = g*(1-(arg(1,1)/g)^0.5-arg(1,1)/g*sin(10*pi*arg(1,1)));
        case 6
            f2 = g*(1-(f1/g)^2);
    end
end


  • 3
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
NSGA-II(非支配排序遗传算法II)是一种经典的多目标优化算法。为了改进NSGA-II算法,在算法的选择、交叉和变异操作中引入了一些改进策略,以提高算法的性能和收敛速度。 在选择操作中,可以采用非支配排序与拥挤度计算相结合的方式来选择个体。非支配排序可以根据个体的非支配等级来判断其优劣,拥挤度计算可以评估个体在目标空间的分布情况,从而保证选择中既考虑到了个体的优越性,又考虑到了多样性。 在交叉操作中,可以采用模拟二进制杂交(SBX)算子,该算子能够在保证交换基因的多样性的同时,保持了较好的搜索效果。通过调整交叉概率和交叉分布指数参数,可以控制交叉操作的强度。 在变异操作中,可以采用多项式变异算子,该算子能够在演化过程中维持一定的多样性,从而获取更多的解集。通过调整变异概率和变异分布指数参数,可以控制变异操作的强度。 为了验证改进的NSGA-II算法,可以使用ZDT(Zitzler-Deb-Thiele)测试函数集进行测试ZDT函数集是一组常用于评估多目标优化算法性能的标准测试函数。通过在Matlab中实现改进后的NSGA-II算法,并将其应用于ZDT函数集,可以对算法的优化效果进行评估和比较。 在实验中,可以通过比较改进NSGA-II算法得到的Pareto前沿集和真实前沿集之间的距离来评估算法的性能。距离越小说明算法的收敛性和准确性越好。同时,还可以比较算法在不同测试函数上的表现,从而进一步分析算法的优劣。 总之,通过改进NSGA-II算法并使用ZDT测试函数集在Matlab中进行验证,可以评估改进算法的性能,并对其进行比较和分析。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值