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