NSGA2 演算法Matlab實現

為了能隨時瞭解Matlab主要操作及思想。

故本文貼上NSGA-Ⅱ演算法Matlab實現（測試函式為ZDT1）。
更多內容訪問omegaxyz.com
NSGA-Ⅱ就是在第一代非支配排序遺傳演算法的基礎上改進而來，其改進主要是針對如上所述的三個方面：
①提出了快速非支配排序演算法，一方面降低了計算的複雜度，另一方面它將父代種群跟子代種群進行合併，使得下一代的種群從雙倍的空間中進行選取，從而保留了最為優秀的所有個體；
②引進精英策略，保證某些優良的種群個體在進化過程中不會被丟棄，從而提高了優化結果的精度；
③採用擁擠度和擁擠度比較運算元，不但克服了NSGA中需要人為指定共享引數的缺陷，而且將其作為種群中個體間的比較標準，使得準Pareto域中的個體能均勻地擴充套件到整個Pareto域，保證了種群的多樣性。

Matlab實現：

function NSGAII()
clc;format compact;tic;hold on

%---初始化/引數設定
    generations=100;                                %迭代次數
    popnum=100;                                     %種群大小(須為偶數)
    poplength=30;                                   %個體長度
    minvalue=repmat(zeros(1,poplength),popnum,1);   %個體最小值
    maxvalue=repmat(ones(1,poplength),popnum,1);    %個體最大值    
    population=rand(popnum,poplength).*(maxvalue-minvalue)+minvalue;    %產生新的初始種群

%---開始迭代進化
    for gene=1:generations                      %開始迭代

%-------交叉 
        newpopulation=zeros(popnum,poplength);  %子代種群
        for i=1:popnum/2                        %交叉產生子代
            k=randperm(popnum);                 %從種群中隨機選出兩個父母,不採用二進位制聯賽方法
            beta=(-1).^round(rand(1,poplength)).*abs(randn(1,poplength))*1.481; %採用正態分佈交叉產生兩個子代
            newpopulation(i*2-1,:)=(population(k(1),:)+population(k(2),:))/2+beta.*(population(k(1),:)-population(k(2),:))./2;    %產生第一個子代           
            newpopulation(i*2,:)=(population(k(1),:)+population(k(2),:))/2-beta.*(population(k(1),:)-population(k(2),:))./2;      %產生第二個子代
        end

%-------變異        
        k=rand(size(newpopulation));    %隨機選定要變異的基因位
        miu=rand(size(newpopulation));  %採用多項式變異
        temp=k<1/poplength & miu<0.5;   %要變異的基因位
        newpopulation(temp)=newpopulation(temp)+(maxvalue(temp)-minvalue(temp)).*((2.*miu(temp)+(1-2.*miu(temp)).*(1-(newpopulation(temp)-minvalue(temp))./(maxvalue(temp)-minvalue(temp))).^21).^(1/21)-1);        %變異情況一
        newpopulation(temp)=newpopulation(temp)+(maxvalue(temp)-minvalue(temp)).*(1-(2.*(1-miu(temp))+2.*(miu(temp)-0.5).*(1-(maxvalue(temp)-newpopulation(temp))./(maxvalue(temp)-minvalue(temp))).^21).^(1/21));  %變異情況二

%-------越界處理/種群合併        
        newpopulation(newpopulation>maxvalue)=maxvalue(newpopulation>maxvalue); %子代越上界處理
        newpopulation(newpopulation<minvalue)=minvalue(newpopulation<minvalue); %子代越下界處理
        newpopulation=[population;newpopulation];                               %合併父子種群

%-------計算目標函式值        
        functionvalue=zeros(size(newpopulation,1),2);           %合併後種群的各目標函式值,這裡的問題是ZDT1
        functionvalue(:,1)=newpopulation(:,1);                  %計算第一維目標函式值
        g=1+9*sum(newpopulation(:,2:poplength),2)./(poplength-1);
        functionvalue(:,2)=g.*(1-(newpopulation(:,1)./g).^0.5); %計算第二維目標函式值

%-------非支配排序        
        fnum=0;                                             %當前分配的前沿面編號
        cz=false(1,size(functionvalue,1));                  %記錄個體是否已被分配編號
        frontvalue=zeros(size(cz));                         %每個個體的前沿面編號
        [functionvalue_sorted,newsite]=sortrows(functionvalue);    %對種群按第一維目標值大小進行排序
        while ~all(cz)                                      %開始迭代判斷每個個體的前沿面,採用改進的deductive sort
            fnum=fnum+1;
            d=cz;
            for i=1:size(functionvalue,1)
                if ~d(i)
                    for j=i+1:size(functionvalue,1)
                        if ~d(j)
                            k=1;                            
                            for m=2:size(functionvalue,2)
                                if functionvalue_sorted(i,m)>functionvalue_sorted(j,m)
                                    k=0;
                                    break
                                end
                            end
                            if k
                                d(j)=true;
                            end
                        end
                    end
                    frontvalue(newsite(i))=fnum;
                    cz(i)=true;
                end
            end
        end

%-------計算擁擠距離/選出下一代個體        
        fnum=0;                                                                 %當前前沿面
        while numel(frontvalue,frontvalue<=fnum+1)<=popnum                      %判斷前多少個面的個體能完全放入下一代種群
            fnum=fnum+1;
        end        
        newnum=numel(frontvalue,frontvalue<=fnum);                              %前fnum個面的個體數
        population(1:newnum,:)=newpopulation(frontvalue<=fnum,:);               %將前fnum個面的個體複製入下一代                       
        popu=find(frontvalue==fnum+1);                                          %popu記錄第fnum+1個面上的個體編號
        distancevalue=zeros(size(popu));                                        %popu各個體的擁擠距離
        fmax=max(functionvalue(popu,:),[],1);                                   %popu每維上的最大值
        fmin=min(functionvalue(popu,:),[],1);                                   %popu每維上的最小值
        for i=1:size(functionvalue,2)                                           %分目標計算每個目標上popu各個體的擁擠距離
            [~,newsite]=sortrows(functionvalue(popu,i));
            distancevalue(newsite(1))=inf;
            distancevalue(newsite(end))=inf;
            for j=2:length(popu)-1
                distancevalue(newsite(j))=distancevalue(newsite(j))+(functionvalue(popu(newsite(j+1)),i)-functionvalue(popu(newsite(j-1)),i))/(fmax(i)-fmin(i));
            end
        end                                      
        popu=-sortrows(-[distancevalue;popu]')';                                %按擁擠距離降序排序第fnum+1個面上的個體
        population(newnum+1:popnum,:)=newpopulation(popu(2,1:popnum-newnum),:); %將第fnum+1個面上擁擠距離較大的前popnum-newnum個個體複製入下一代        
    end

%---程式輸出    
    fprintf('已完成,耗時%4s秒\n',num2str(toc));          %程式最終耗時
    output=sortrows(functionvalue(frontvalue==1,:));    %最終結果:種群中非支配解的函式值
    plot(output(:,1),output(:,2),'*b');                 %作圖
    axis([0,1,0,1]);xlabel('F_1');ylabel('F_2');title('ZDT1')
end

更多內容訪問omegaxyz.com

網站文章採用知識共享許可協議BY-NC-SA4.0授權
© 2018 • OmegaXYZ-版權所有 轉載請註明出處