nsga-II算法的matlab代码

以下是NSGA-II算法的MATLAB代码：

%% NSGA-II Algorithm
function [pop, front] = NSGAII(pop_size,gen_size,vari_num,x_min,x_max)

% Parameters Setting
pc = 0.9;   % crossover probability
pm = 1/vari_num; % mutation probability
eta_c = 20; % crossover distribution index
eta_m = 20; % mutation distribution index
pop = rand(pop_size,vari_num) .* repmat((x_max-x_min),pop_size,1) + repmat(x_min,pop_size,1); % initial population
for i = 1:pop_size
    obj(i,1) = f1(pop(i,:));
    obj(i,2) = f2(pop(i,:));
end

% NSGA-II Algorithm
for i = 1:gen_size
    % Non-dominated Sorting
    [front,~] = non_domination_sort(obj);
    % Crowding Distance Calculation
    for j = 1:length(front)
        crowd_dis(j,:) = crowding_distance(obj(front{j},:));
    end
    % Mating Selection & Variation
    pop_new = zeros(pop_size,vari_num);
    obj_new = zeros(pop_size,2);
    count = 0;
    for j = 1:length(front)
        [temp,index] = sort(crowd_dis(j,:),'descend');
        front_member = front{j};
        pop_temp = pop(front_member,:);
        for k = 1:length(front_member)
            if rand < pc && k ~= length(front_member)
                p1 = front_member(k);
                p2 = front_member(k+1);
                pop_new(count+1,:) = crossover(pop_temp(p1,:),pop_temp(p2,:),eta_c);
                obj_new(count+1,:) = [f1(pop_new(count+1,:)),f2(pop_new(count+1,:))];
                count = count + 1;
            elseif rand < pm
                p = front_member(k);
                pop_new(count+1,:) = mutation(pop_temp(p,:),eta_m,x_min,x_max);
                obj_new(count+1,:) = [f1(pop_new(count+1,:)),f2(pop_new(count+1,:))];
                count = count + 1;
            end
        end
        if count >= pop_size
            break;
        end
    end
    % Combine Parent & Offspring Populations
    pop = [pop;pop_new(1:pop_size-count,:)];
    obj = [obj;obj_new(1:pop_size-count,:)];
end

% Results Output
[front,~] = non_domination_sort(obj);
for i = 1:length(front)
    plot(obj(front{i},1),obj(front{i},2),'o');
    hold on;
end
xlabel('f_1');
ylabel('f_2');
end

%% Non-dominated Sorting
function [front,rank] = non_domination_sort(obj)

n = size(obj,1);
rank = inf(1,n);
dominate = false(n);

for i = 1:n
    for j = i+1:n
        if all(obj(i,:) <= obj(j,:)) && any(obj(i,:) < obj(j,:))
            dominate(i,j) = true;
        elseif all(obj(i,:) >= obj(j,:)) && any(obj(i,:) > obj(j,:))
            dominate(j,i) = true;
        end
    end
end

for i = 1:n
    if sum(dominate(:,i)) == 0
        rank(i) = 1;
    end
end

front{1} = find(rank == 1);

i = 1;
while ~isempty(front{i})
    Q = [];
    for j = front{i}
        for k = 1:n
            if dominate(j,k)
                dominate(j,k) = false;
                if sum(dominate(k,:)) == 0
                    rank(k) = i + 1;
                    Q = [Q,k];
                end
            end
        end
    end
    i = i + 1;
    front{i} = Q;
end
end

%% Crowding Distance Calculation
function [crowd_dis] = crowding_distance(obj)

n = size(obj,1);
crowd_dis = zeros(1,n);
f_max = max(obj,[],1);
f_min = min(obj,[],1);

for i = 1:size(obj,2)
    [~,index] = sort(obj(:,i));
    crowd_dis(index(1)) = inf;
    crowd_dis(index(end)) = inf;
    for j = 2:n-1
        crowd_dis(index(j)) = crowd_dis(index(j)) + (obj(index(j+1),i) - obj(index(j-1),i))/(f_max(i) - f_min(i));
    end
end
end

%% SBX Crossover Operator
function [offspring] = crossover(p1,p2,eta_c)

n = length(p1);
u = rand(1,n);
offspring = zeros(1,n);

for i = 1:n
    if u(i) <= 0.5
        if abs(p1(i)-p2(i)) > 1e-10
            if p1(i) < p2(i)
                x1 = p1(i);
                x2 = p2(i);
            else
                x1 = p2(i);
                x2 = p1(i);
            end
            y1 = (x1 - floor(x1)) + floor(x2);
            y2 = (x2 - floor(x2)) + floor(x1);
            if rand < 0.5
                offspring(i) = y1;
            else
                offspring(i) = y2;
            end
        else
            offspring(i) = p1(i);
        end
    else
        if abs(p1(i)-p2(i)) > 1e-10
            if p1(i) < p2(i)
                x1 = p1(i);
                x2 = p2(i);
            else
                x1 = p2(i);
                x2 = p1(i);
            end
            y1 = 2*x1 - x2;
            y2 = 2*x2 - x1;
            if y1 < 0
                y1 = 0;
            elseif y1 > 1
                y1 = 1;
            end
            if y2 < 0
                y2 = 0;
            elseif y2 > 1
                y2 = 1;
            end
            offspring(i) = y1;
        else
            offspring(i) = p1(i);
        end
    end
end
end

%% Polynomial Mutation Operator
function [offspring] = mutation(p,eta_m,x_min,x_max)

n = length(p);
offspring = p;

for i = 1:n
    if rand < eta_m/n
        delta1 = (p(i) - x_min)/(x_max - x_min);
        delta2 = (x_max - p(i))/(x_max - x_min);
        u = rand;
        if u <= 0.5
            deltaq = (2*u + (1 - 2*u)*(1-delta1)^(eta_m+1))^(1/(eta_m+1))-1;
        else
            deltaq = 1 - (2*(1-u) + 2*(u-0.5)*(1-delta2)^(eta_m+1))^(1/(eta_m+1));
        end
        offspring(i) = p(i) + deltaq*(x_max - x_min);
        if offspring(i) < x_min
            offspring(i) = x_min;
        elseif offspring(i) > x_max
            offspring(i) = x_max;
        end
    end
end
end

%% Objective Function 1
function [y1] = f1(x)

y1 = x(1);

end

%% Objective Function 2
function [y2] = f2(x)

y2 = (1+x(2))/x(1);

end

其中，`pop_size`为种群大小，`gen_size`为迭代次数，`vari_num`为变量个数，`x_min`和`x_max`分别为变量的上下界。`f1`和`f2`分别为两个目标函数，根据具体问题进行修改。