基于非支配排序的多目标PSO算法MATLAB实现

function [solution,time,pop,pfitness,site,LeaderAVE] = NSPSO(train_F,train_L)
tic
global maxFES
dim = size(train_F,2);
FES = 1;
sizep = 30;
pop = rand(sizep,dim);
popv = rand(sizep,dim);
pfitness = zeros(sizep,2);
LeaderAVE = zeros(1,2);
while FES <maxFES
    Off_P = zeros(sizep,dim);
    Off_V = zeros(sizep,dim);
    ofitness = zeros(sizep,2);
    for i=1:sizep
        [pfitness(i,1),pfitness(i,2)] = FSKNN(pop(i,:),i,train_F,train_L);
    end
    Front = NDSort(pfitness(:,1:2),sizep);
    [~,rank] = sortrows([Front',-CrowdingDistance(pfitness,Front)']);
    LeaderSet = rank(1:10);
    solution = pfitness(LeaderSet,:);
    LeaderAVE(1) = mean(solution(:,1));
    LeaderAVE(2) = mean(solution(:,2));
    for i = 1:sizep
        good = LeaderSet(randperm(length(LeaderSet),1));
        r1 = rand(1,dim);
        r2 = rand(1,dim);
        Off_V(i,:) = r1.*popv(i,:) +  r2.*(pop(good,:)-pop(i,:));
        Off_P(i,:) = pop(i,:) + Off_V(i,:);
    end
    for i=1:sizep
    	[ofitness(i,1),ofitness(i,2)] = FSKNN(Off_P(i,:),i,train_F,train_L);
    end
    temppop = [pop;Off_P];
    tempv = [popv;Off_V];
    tempfiness = [pfitness;ofitness];
    [FrontNO,MaxFNO] = NDSort(tempfiness(:,1:2),sizep);
    Next = false(1,length(FrontNO));
    Next(FrontNO<MaxFNO) = true;  
    PopObj = tempfiness;
    fmax   = max(PopObj(FrontNO==1,:),[],1);
    fmin   = min(PopObj(FrontNO==1,:),[],1);
    PopObj = (PopObj-repmat(fmin,size(PopObj,1),1))./repmat(fmax-fmin,size(PopObj,1),1);
    
    % Select the solutions in the last front
    Last = find(FrontNO==MaxFNO);
    del  = Truncation(PopObj(Last,:),length(Last)-sizep+sum(Next));
    Next(Last(~del)) = true;
    % Population for next generation
    pop = temppop(Next,:);
    popv = tempv(Next,:);
    pfitness = tempfiness(Next,:);
    fprintf('GEN: %2d   Error: %.4f  F:%.2f\n',FES,LeaderAVE(1),LeaderAVE(2));
    FES = FES + 1;
end
[FrontNO,~] = NDSort(pfitness(:,1:2),sizep);
site = find(FrontNO==1);
solution = pfitness(site,:);
LeaderAVE(1) = mean(solution(:,1));
LeaderAVE(2) = mean(solution(:,2));
toc
time = toc;
end

function CrowdDis = CrowdingDistance(PopObj,FrontNO)
% Calculate the crowding distance of each solution front by front

% Copyright 2015-2016 Ye Tian

    [N,M]    = size(PopObj);
    CrowdDis = zeros(1,N);
    Fronts   = setdiff(unique(FrontNO),inf);
    for f = 1 : length(Fronts)
        Front = find(FrontNO==Fronts(f));
        Fmax  = max(PopObj(Front,:),[],1);
        Fmin  = min(PopObj(Front,:),[],1);
        for i = 1 : M
            [~,Rank] = sortrows(PopObj(Front,i));
            CrowdDis(Front(Rank(1)))   = inf;
            CrowdDis(Front(Rank(end))) = inf;
            for j = 2 : length(Front)-1
                CrowdDis(Front(Rank(j))) = CrowdDis(Front(Rank(j)))+(PopObj(Front(Rank(j+1)),i)-PopObj(Front(Rank(j-1)),i))/(Fmax(i)-Fmin(i));
            end
        end
    end
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    N = size(PopObj,1);
    
    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,N);
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end