多目标遗传算法SPEA2的matlab代码

以下是多目标遗传算法SPEA2的Matlab代码：

% SPEA2算法
function [pop,popObj] = SPEA2(Global, pop)
% 参数设置
N = Global.N;
K = ceil(sqrt(N));
T = 20;
maxgen = Global.maxgen;
% 计算适应度值
popObj = pop.objs;
[n,~] = size(popObj);
D = pdist2(popObj,popObj);
D(D==0) = inf;
[~,knn] = sort(D,2);
knn = knn(:,1:K);
S = zeros(n,K);
for i = 1 : n
    S(i,:) = unique([knn(i,:),i]);
end
fit = zeros(1,n);
for i = 1 : n
    fit(i) = 1 / (mean(min(D(S(i,:),S(i,:)),[],2)) + 1e-6);
end
% 计算拥挤度
dis = zeros(1,n);
for i = 1 : Global.M
    [~,rank] = sort(popObj(:,i));
    dis(rank(1))   = inf;
    dis(rank(end)) = inf;
    for j = 2 : n-1
        dis(rank(j)) = dis(rank(j)) + (popObj(rank(j+1),i)-popObj(rank(j-1),i))/(popObj(rank(end),i)-popObj(rank(1),i)+eps);
    end
end
% 进化
gen = 1;
while gen <= maxgen
    MatingPool = TournamentSelection(2, N, fit);
    Offspring  = GeneticVariation(pop(MatingPool));
    popAll     = [pop, Offspring];
    popObjAll  = [popAll.objs];
    [nAll,~]   = size(popObjAll);
    D = pdist2(popObjAll,popObjAll);
    D(D==0) = inf;
    [~ ,knn] = sort(D,2);
    knn = knn(:,1:K);
    S = zeros(nAll,K);
    for i = 1:nAll
        S(i,:) = unique([knn(i,:),i]);
    end
    fitAll = zeros(1,nAll);
    for i = 1:nAll
        fitAll(i) = 1 / (mean(min(D(S(i,:),S(i,:)),[],2)) + 1e-6);
    end
    objSort = sortrows([popObjAll;popObj]);
    disSort = [dis,dis];
    N = min(T,length(fitAll));
    Next = EnvironmentalSelection([popAll, pop], [fitAll, fit], objSort(1:N,:), disSort(1:N));
    % 更新种群
    pop    = Next(1:Global.N);
    popObj = pop.objs;
%     fprintf('gen=%d, f1=%e\n',gen,min(popObj(:,1)))
    gen = gen + 1;
end
end

function Population = TournamentSelection(N, T, Fitness)
% Tournament selection
% N: the number of selected solutions
% T: the tournament size
% Fitness: the fitness values of the current population
% Population: the indices of the selected solutions
Population = zeros(N,1);
for i = 1 : N
    Tournament = randperm(T,T);
    [~,best] = min(Fitness(Tournament));
    Population(i) = Tournament(best);
end
end

function Offspring = GeneticVariation(Parent)
% Genetic variation operators
% Parent: the current population
% Offspring: the offspring population
    [N,D] = size(Parent(1).vars);
    Parent = Parent(randperm(length(Parent)));
    Offspring = Parent;
    for i = 1 : 2 : length(Parent)
        p1 = Parent(i);
        p2 = Parent(i+1);
        if rand < 0.9 % crossover
            % Simulated binary crossover
            beta = zeros(1,D);
            mu   = rand(1,D) < 0.5;
            beta(mu) = (2*rand(1,sum(mu))).^(1/(Global.V+1));
            beta(~mu)= (2*rand(1,sum(~mu))).^(-1/(Global.V+1));
            beta = beta.*(-1).^randi([0,1],1,D);
            c1 = 0.5*((1+beta).*p1.vars + (1-beta).*p2.vars);
            c2 = 0.5*((1-beta).*p1.vars + (1+beta).*p2.vars);
        else % mutation
            % polynomial mutation
            Site = rand(N,D) < 1/D;
            mu   = rand(N,D);
            temp = Site & mu<=0.5;
            Offspring(i).vars(temp) = p1.vars(temp)+(p1.vars(temp)-p2.vars(temp)).*(2.*mu(temp)+(1-2.*mu(temp)).*...
                (abs((p1.vars(temp)-p2.vars(temp))./(1e-16+1-p1.vars(temp)+p2.vars(temp)))));
            temp = Site & mu>0.5;
            Offspring(i).vars(temp) = p2.vars(temp)+(p2.vars(temp)-p1.vars(temp)).*(2.*(mu(temp)-0.5)+(1-2.*(mu(temp)-0.5)).*...
                (abs((p1.vars(temp)-p2.vars(temp))./(1e-16+1-p1.vars(temp)+p2.vars(temp)))));
            temp = ~Site;
            Offspring(i).vars(temp) = p1.vars(temp);
            c1 = Offspring(i).vars;
            Site = rand(1,D) < 1/D;
            mu   = rand(1,D);
            temp = Site & mu<=0.5;
            Offspring(i+1).vars(temp) = p2.vars(temp)+(p2.vars(temp)-p1.vars(temp)).*(2.*mu(temp)+(1-2.*mu(temp)).*...
                (abs((p1.vars(temp)-p2.vars(temp))./(1e-16+1-p1.vars(temp)+p2.vars(temp)))));
            temp = Site & mu>0.5;
            Offspring(i+1).vars(temp) = p1.vars(temp)+(p1.vars(temp)-p2.vars(temp)).*(2.*(mu(temp)-0.5)+(1-2.*(mu(temp)-0.5)).*...
                (abs((p1.vars(temp)-p2.vars(temp))./(1e-16+1-p1.vars(temp)+p2.vars(temp)))));
            temp = ~Site;
            Offspring(i+1).vars(temp) = p2.vars(temp);
            c2 = Offspring(i+1).vars;
        end
        Offspring(i).obj = Global.evaluate(c1);
        Offspring(i+1).obj = Global.evaluate(c2);
    end
end

function Population = EnvironmentalSelection(Population, Fitness, Obj, Dis)
% Environmental selection
% Population: the current population
% Fitness: the fitness values of the current population
% Obj: the objective values of the whole population
% Dis: the crowding distance of the whole population
% Population: the population after environmental selection
    N = length(Population);
    M = size(Obj,2);
    [~,Rank] = sort(Fitness);
    Population = Population(Rank);
    Obj = Obj(Rank,:);
    Dis = Dis(Rank);
    Next = [];
    LastFront = 1;
    while length(Next) + length(Obj(LastFront:end,1)) <= N
        Next = [Next, Population(LastFront : end)];
        LastFront = find(cumsum(Dis(LastFront:end))/sum(Dis(LastFront:end))>=1,1) + LastFront - 1;
    end
    if length(Next) < N
        [~,Rank] = sort(Obj(LastFront:end,:),1);
        Rank = Rank + LastFront - 1;
        Front = unique(Rank(:,1));
        Remain = N - length(Next);
        for i = 1 : length(Front)
            if Remain == 0
                break;
            end
            FrontObj = Rank(Rank(:,1)==Front(i),2:end);
            [~,Rank2] = sort(Dis(FrontObj),'descend');
            Next = [Next, Population(FrontObj(Rank2(1:min(length(Rank2),Remain))))];
            Remain = N - length(Next);
        end
    end
end

这段代码可以配合Matlab上的优化工具箱使用，通过`Global.evaluate`函数来计算适应度值。需要注意的是，`Global`变量是一个结构体，需要在调用函数时进行初始化。