多目標優化演算法（一）NSGA-Ⅱ（NSGA2）

曉風wangchao發表於2018-09-28
優化演算法

多目標優化演算法（一）NSGA-Ⅱ(NSGA2)

注：沒有想到這篇部落格竟然有很多人檢視，這是我去年寫的演算法，裡面難免會有一些錯誤，大家可以看看評論區，這裡的matlab程式碼寫的不太好，是以C語言思維寫的，基本上沒有並行，初學者建議可以看看platEMO上的原始碼，提前培養好寫程式碼的習慣！

0. 前言

這個演算法是本人接觸科研學習實現的第一個演算法，因此想在這裡和大家分享一下心得。

1. 演算法簡介

NSGA-Ⅱ演算法，即帶有精英保留策略的快速非支配多目標優化演算法，是一種基於Pareto最優解的多目標優化演算法。

1.1 Pareto支配關係以及Pareto等級

Pareto支配關係：對於最小化多目標優化問題，對於n個目標分量 fi(x),i=1...nf_i(x), i=1...n，任意給定兩個決策變數XaX_a，XbX_b,如果有以下兩個條件成立，則稱XaX_a支配XbX_b。
1.對於∀i∈1,2,...,n\forall i \in {1,2,...,n}，都有fi(Xa)≤fi(Xb)f_i(X_a) \leq f_i(X_b)成立。
2.∃i∈1,2,...,n\exists i \in {1,2,...,n}，使得fi(Xa)<fi(Xb)f_i(X_a) <f_i(X_b)成立。
如果對於一個決策變數，不存在其他決策變數能夠支配他，那麼就稱該決策變數為非支配解。
Pareto等級：在一組解中，非支配解Pareto等級定義為1，將非支配解從解的集合中刪除，剩下解的Pareto等級定義為2，依次類推，可以得到該解集合中所有解的Pareto等級。示意圖如圖1所示。
在這裡插入圖片描述

1.2 快速非支配排序

假設種群大小為P，該演算法需要計算每個個體p的被支配個數npn_p和該個體支配的解的集合SpS_p這兩個引數。遍歷整個種群，該引數的計算複雜度為O(mN2)O(mN^2)。該演算法的虛擬碼如下：
1.計算出種群中每個個體的兩個引數npn_p和SpS_p。
2.將種群中引數np=0n_p=0的個體放入集合F1F_1中。
3.for 個體i∈F1i \in F_1:
        for 個體l∈Sil \in S_i：
               nl=nl−1n_l=n_l-1；(該步驟即消除Pareto等級1對其餘個體的支配，相當於將Pareto等級1的個體從集合中刪除)
                if nl=0n_l=0：
                       將個體l加入集合F2F_2中。
                end
        end
end
4.上面得到Pareto等級2的個體的集合F2F_2，對集合F2F_2中的個體繼續重複步驟3，依次類推直到種群等級被全部劃分。
matlab 程式碼如下：

function [F,chromo] = non_domination_sort( pop,chromo,f_num,x_num )
%non_domination_sort 初始種群的非支配排序和計算擁擠度
%初始化pareto等級為1
pareto_rank=1;
F(pareto_rank).ss=[];%pareto等級為pareto_rank的集合
p=[];%每個個體p的集合
for i=1:pop
    %%%計算出種群中每個個體p的被支配個數n和該個體支配的解的集合s
    p(i).n=0;%被支配個體數目n
    p(i).s=[];%支配解的集合s
    for j=1:pop
        less=0;%y'的目標函式值小於個體的目標函式值數目
        equal=0;%y'的目標函式值等於個體的目標函式值數目
        greater=0;%y'的目標函式值大於個體的目標函式值數目
        for k=1:f_num
            if(chromo(i,x_num+k)<chromo(j,x_num+k))
                less=less+1;
            elseif(chromo(i,x_num+k)==chromo(j,x_num+k))
                equal=equal+1;
            else
                greater=greater+1;
            end
        end
        if(less==0 && equal~=f_num)%if(less==0 && greater~=0)
            p(i).n=p(i).n+1;%被支配個體數目n+1
        elseif(greater==0 && equal~=f_num)%elseif(greater==0 && less~=0)
            p(i).s=[p(i).s j];
        end
    end
    %%%將種群中引數為n的個體放入集合F(1)中
    if(p(i).n==0)
        chromo(i,f_num+x_num+1)=1;%儲存個體的等級資訊
        F(pareto_rank).ss=[F(pareto_rank).ss i];
    end
end
%%%求pareto等級為pareto_rank+1的個體
while ~isempty(F(pareto_rank).ss)
    temp=[];
    for i=1:length(F(pareto_rank).ss)
        if ~isempty(p(F(pareto_rank).ss(i)).s)
            for j=1:length(p(F(pareto_rank).ss(i)).s)
                p(p(F(pareto_rank).ss(i)).s(j)).n=p(p(F(pareto_rank).ss(i)).s(j)).n - 1;%nl=nl-1
                if p(p(F(pareto_rank).ss(i)).s(j)).n==0
                    chromo(p(F(pareto_rank).ss(i)).s(j),f_num+x_num+1)=pareto_rank+1;%儲存個體的等級資訊
                    temp=[temp p(F(pareto_rank).ss(i)).s(j)];
                end
            end
        end
    end
    pareto_rank=pareto_rank+1;
    F(pareto_rank).ss=temp;
end

1.3 擁擠度

為了使得到的解在目標空間中更加均勻，這裡引入了擁擠度ndn_d，其演算法如下：

  1. 令引數nd=0，n∈1…Nn_d=0，n∈1…N。
  2. for 每個目標函式fmf_m ：
            ① 根據該目標函式對該等級的個體進行排序，記fmmaxf_m^{max}為個體目標函式值fmf_m的最大值，fmminf_m^{min}為個體目標函式值fmf_m的最小值；
            ② 對於排序後兩個邊界的擁擠度1d1_d和NdN_d置為∞；
            ③ 計算nd=nd+(fm(i+1)−fm(i−1))/(fmmax−fmmin)n_d=n_d+(f_m (i+1)-f_m (i-1))/(f_m^{max}-f_m^{min}),其中fm(i+1)f_m (i+1)是該個體排序後後一位的目標函式值。

    end

      從二目標優化問題來看，就像是該個體在目標空間所能生成的最大的矩形（該矩形不能觸碰目標空間其他的點）的邊長之和。擁擠度示意圖如圖2所示。
擁擠度示意圖
matlab程式碼如下：

function chromo = crowding_distance_sort( F,chromo,f_num,x_num )
%計算擁擠度
%%%按照pareto等級對種群中的個體進行排序
[~,index]=sort(chromo(:,f_num+x_num+1));
[~,mm1]=size(chromo);
temp=zeros(length(index),mm1);
for i=1:length(index)%=pop
    temp(i,:)=chromo(index(i),:);%按照pareto等級排序後種群
end

%%%對於每個等級的個體開始計算擁擠度
current_index = 0;
for pareto_rank=1:(length(F)-1)%計算F的迴圈時多了一次空，所以減掉
    %%擁擠度初始化為0
    nd=[];
    nd(:,1)=zeros(length(F(pareto_rank).ss),1);
    %y=[];%儲存當前處理的等級的個體
    [~,mm2]=size(temp);
    y=zeros(length(F(pareto_rank).ss),mm2);%儲存當前處理的等級的個體
    previous_index=current_index + 1;
    for i=1:length(F(pareto_rank).ss)
        y(i,:)=temp(current_index + i,:);
    end
    current_index=current_index + i;
    %%對於每一個目標函式fm
    for i=1:f_num
        %%根據該目標函式值對該等級的個體進行排序
        [~,index_objective]=sort(y(:,x_num+i));
        %objective_sort=[];%通過目標函式排序後的個體
        [~,mm3]=size(y);
        objective_sort=zeros(length(index_objective),mm3);%通過目標函式排序後的個體
        for j=1:length(index_objective)
            objective_sort(j,:)=y(index_objective(j),:);
        end
        %%記fmax為最大值，fmin為最小值
        fmin=objective_sort(1,x_num+i);
        fmax=objective_sort(length(index_objective),x_num+i);
        %%對排序後的兩個邊界擁擠度設為1d和nd設為無窮
        y(index_objective(1),f_num+x_num+1+i)=Inf;
        y(index_objective(length(index_objective)),f_num+x_num+1+i)=Inf;
        %%計算nd=nd+(fm(i+1)-fm(i-1))/(fmax-fmin)
        for j=2:(length(index_objective)-1)
            pre_f=objective_sort(j-1,x_num+i);
            next_f=objective_sort(j+1,x_num+i);
            if (fmax-fmin==0)
                y(index_objective(j),f_num+x_num+1+i)=Inf;
            else
                y(index_objective(j),f_num+x_num+1+i)=(next_f-pre_f)/(fmax-fmin);
            end
        end
    end
    %多個目標函式擁擠度求和
    for i=1:f_num
        nd(:,1)=nd(:,1)+y(:,f_num+x_num+1+i);
    end
    %第2列儲存擁擠度，其他的覆蓋掉
    y(:,f_num+x_num+2)=nd;
    y=y(:,1:(f_num+x_num+2));
    temp_two(previous_index:current_index,:) = y;
end
chromo=temp_two;
end

1.4 精英保留策略

1.首先將父代種群CiC_i和子代種群DiD_i合成種群RiR_i。
2.根據以下規則從種群RiR_i生成新的父代種群Ci+1C_{i+1}：
      ①根據Pareto等級從低到高的順序，將整層種群放入父代種群Ci+1C_{i+1}，直到某一層該層個體不能全部放入父代種群Ci+1C_{i+1}；
      ②將該層個體根據擁擠度從大到小排列，依次放入父代種群Ci+1C_{i+1}中，直到父代種群Ci+1C_{i+1}填滿。
matlab程式碼如下：

function chromo = elitism( pop,combine_chromo2,f_num,x_num )
%精英保留策略
[pop2,temp]=size(combine_chromo2);
chromo_rank=zeros(pop2,temp);
chromo=zeros(pop,(f_num+x_num+2));
%根據pareto等級從高到低進行排序
[~,index_rank]=sort(combine_chromo2(:,f_num+x_num+1));
for i=1:pop2
    chromo_rank(i,:)=combine_chromo2(index_rank(i),:);
end
%求出最高的pareto等級
max_rank=max(combine_chromo2(:,f_num+x_num+1));
%根據排序後的順序，將等級相同的種群整個放入父代種群中，直到某一層不能全部放下為止
prev_index=0;
for i=1:max_rank
    %尋找當前等級i個體裡的最大索引
    current_index=max(find(chromo_rank(:,(f_num+x_num+1))==i));
    %不能放下為止
    if(current_index>pop)
        %剩餘群體數
        remain_pop=pop-prev_index;
        %等級為i的個體
        temp=chromo_rank(((prev_index+1):current_index),:);
        %根據擁擠度從大到小排序
        [~,index_crowd]=sort(temp(:,(f_num+x_num+2)),'descend');
        %填滿父代
        for j=1:remain_pop
            chromo(prev_index+j,:)=temp(index_crowd(j),:);
        end
        return;
    elseif(current_index<pop)
        % 將所有等級為i的個體直接放入父代種群
        chromo(((prev_index+1):current_index),:)=chromo_rank(((prev_index+1):current_index),:);
    else
        % 將所有等級為i的個體直接放入父代種群
        chromo(((prev_index+1):current_index),:)=chromo_rank(((prev_index+1):current_index),:);
        return;
    end
    prev_index = current_index;
end
end

1.5 實數編碼的交叉操作（SBX）

模擬二進位制交叉：
x1j(t)=0.5×[(1+γj)x1j(t)+(1−γj)x2j(t)]x _{1j}(t)=0.5×[(1+γ_j ) x_{1j}(t)+(1-γ_j ) x_{2j} (t)]
x2j(t)=0.5×[(1−γj)x1j(t)+(1+γj)x2j(t)]x _{2j} (t)=0.5×[(1-γ_j ) x_{1j}(t)+(1+γ_j ) x_{2j}(t)]
其中
γj={(2uj)1/(η+1)                        uj<0.5(1/(2(1−uj))1/(η+1)        elseγ_j=\left\{\begin{matrix}(2u_j)^{1/(η+1)}\;\;\;\;\;\;\;\;\;\;\;\;u_j<0.5 \\ (1/(2(1-u_j))^{1/(η+1)}\;\;\;\;else \end{matrix}\right.

1.6 多項式變異（polynomial mutation）

多項式變異：
x1j(t)=x1j(t)+∆jx _{1j} (t)=x_{1j} (t)+∆_j
其中
∆j={(2uj)1/(η+1)−1                        uj<0.5(1−(2(1−uj))1/(η+1)        else∆_j=\left\{\begin{matrix}(2u_j)^{1/(η+1)}-1\;\;\;\;\;\;\;\;\;\;\;\;u_j<0.5 \\ (1-(2(1-u_j))^{1/(η+1)}\;\;\;\;else \end{matrix}\right.
並且0≤uju_j≤1。
matlab程式碼如下：

function chromo_offspring = cross_mutation( chromo_parent,f_num,x_num,x_min,x_max,pc,pm,yita1,yita2,fun )
%模擬二進位制交叉與多項式變異
[pop,~]=size(chromo_parent);
suoyin=1;
for i=1:pop
   %%%模擬二進位制交叉
   %初始化子代種群
   %隨機選取兩個父代個體
   parent_1=round(pop*rand(1));
   if (parent_1<1)
       parent_1=1;
   end
   parent_2=round(pop*rand(1));
   if (parent_2<1)
       parent_2=1;
   end
   %確定兩個父代個體不是同一個
   while isequal(chromo_parent(parent_1,:),chromo_parent(parent_2,:))
       parent_2=round(pop*rand(1));
       if(parent_2<1)
           parent_2=1;
       end
   end
   chromo_parent_1=chromo_parent(parent_1,:);
   chromo_parent_2=chromo_parent(parent_2,:);
   off_1=chromo_parent_1;
   off_2=chromo_parent_1;
   if(rand(1)<pc)
       %進行模擬二進位制交叉
       u1=zeros(1,x_num);
       gama=zeros(1,x_num);
       for j=1:x_num
           u1(j)=rand(1);
           if u1(j)<0.5
               gama(j)=(2*u1(j))^(1/(yita1+1));
           else
               gama(j)=(1/(2*(1-u1(j))))^(1/(yita1+1));
           end
           off_1(j)=0.5*((1+gama(j))*chromo_parent_1(j)+(1-gama(j))*chromo_parent_2(j));
           off_2(j)=0.5*((1-gama(j))*chromo_parent_1(j)+(1+gama(j))*chromo_parent_2(j));
           %使子代在定義域內
           if(off_1(j)>x_max(j))
               off_1(j)=x_max(j);
           elseif(off_1(j)<x_min(j))
               off_1(j)=x_min(j);
           end
           if(off_2(j)>x_max(j))
               off_2(j)=x_max(j);
           elseif(off_2(j)<x_min(j))
               off_2(j)=x_min(j);
           end
       end
       %計運算元代個體的目標函式值
       off_1(1,(x_num+1):(x_num+f_num))=object_fun(off_1,f_num,x_num,fun);
       off_2(1,(x_num+1):(x_num+f_num))=object_fun(off_2,f_num,x_num,fun);
   end
   %%%多項式變異
   if(rand(1)<pm)
       u2=zeros(1,x_num);
       delta=zeros(1,x_num);
       for j=1:x_num
           u2(j)=rand(1);
           if(u2(j)<0.5)
               delta(j)=(2*u2(j))^(1/(yita2+1))-1;
           else
               delta(j)=1-(2*(1-u2(j)))^(1/(yita2+1));
           end
           off_1(j)=off_1(j)+delta(j);
           %使子代在定義域內
           if(off_1(j)>x_max(j))
               off_1(j)=x_max(j);
           elseif(off_1(j)<x_min(j))
               off_1(j)=x_min(j);
           end
       end
       %計運算元代個體的目標函式值
       off_1(1,(x_num+1):(x_num+f_num))=object_fun(off_1,f_num,x_num,fun);
   end
   if(rand(1)<pm)
       u2=zeros(1,x_num);
       delta=zeros(1,x_num);
       for j=1:x_num
           u2(j)=rand(1);
           if(u2(j)<0.5)
               delta(j)=(2*u2(j))^(1/(yita2+1))-1;
           else
               delta(j)=1-(2*(1-u2(j)))^(1/(yita2+1));
           end
           off_2(j)=off_2(j)+delta(j);
           %使子代在定義域內
           if(off_2(j)>x_max(j))
               off_2(j)=x_max(j);
           elseif(off_2(j)<x_min(j))
               off_2(j)=x_min(j);
           end
       end
       %計運算元代個體的目標函式值
       off_2(1,(x_num+1):(x_num+f_num))=object_fun(off_2,f_num,x_num,fun);
   end
   off(suoyin,:)=off_1;
   off(suoyin+1,:)=off_2;
   suoyin=suoyin+2;
end
chromo_offspring=off;
end

1.7 競標賽選擇（tournament selection）

錦標賽法是選擇操作的一種常用方法，二進位制競標賽用的最多。
假設種群規模為N，該法的步驟為：
1.從這N個個體中隨機選擇k（k<n）個個體，k的取值小，效率就高（節省執行時間），但不宜太小，一般取為n/2（取整）。（二進位制即取2）
2.根據每個個體的適應度值，選擇其中適應度值最好的個體進入下一代種群。
3.重複1-2步，至得到新的N個個體。
二進位制選擇的matlab程式碼如下：

function chromo_parent = tournament_selection( chromo )
%二進位制競標賽
[pop, suoyin] = size(chromo);
touranment=2;
a=round(pop/2);
chromo_candidate=zeros(touranment,1);
chromo_rank=zeros(touranment,1);
chromo_distance=zeros(touranment,1);
chromo_parent=zeros(a,suoyin);
% 獲取等級的索引
rank = suoyin - 1;
% 獲得擁擠度的索引
distance = suoyin;
for i=1:a
  for j=1:touranment
      chromo_candidate(j)=round(pop*rand(1));%隨機產生候選者
      if(chromo_candidate(j)==0)%索引從1開始
          chromo_candidate(j)=1;
      end
      if(j>1)
          while (~isempty(find(chromo_candidate(1:j-1)==chromo_candidate(j))))
              chromo_candidate(j)=round(pop*rand(1));
              if(chromo_candidate(j)==0)%索引從1開始
                  chromo_candidate(j)=1;
              end
          end
      end
  end
  for j=1:touranment
      chromo_rank(j)=chromo(chromo_candidate(j),rank);
      chromo_distance(j)=chromo(chromo_candidate(j),distance);
  end
  %取出低等級的個體索引
  minchromo_candidate=find(chromo_rank==min(chromo_rank));
  %多個索引按擁擠度排序
  if (length(minchromo_candidate)~=1)
      maxchromo_candidate=find(chromo_distance(minchromo_candidate)==max(chromo_distance(minchromo_candidate)));
      if(length(maxchromo_candidate)~=1)
          maxchromo_candidate = maxchromo_candidate(1);
      end
      chromo_parent(i,:)=chromo(chromo_candidate(minchromo_candidate(maxchromo_candidate)),:);
  else
      chromo_parent(i,:)=chromo(chromo_candidate(minchromo_candidate(1)),:);
  end
end
end

2. 演算法實現框圖

NSGA-Ⅱ演算法的基本思想為：首先，隨機產生規模為N的初始種群，非支配排序後通過遺傳演算法的選擇、交叉、變異三個基本操作得到第一代子代種群；其次，從第二代開始，將父代種群與子代種群合併，進行快速非支配排序，同時對每個非支配層中的個體進行擁擠度計算，根據非支配關係以及個體的擁擠度選取合適的個體組成新的父代種群；最後，通過遺傳演算法的基本操作產生新的子代種群：依此類推，直到滿足程式結束的條件。該演算法的流程圖如圖3所示。
NSGA-Ⅱ演算法流程圖
初始化matlab程式碼如下：

function chromo = initialize( pop,f_num,x_num,x_min,x_max,fun )
%   種群初始化
for i=1:pop
   for j=1:x_num
       chromo(i,j)=x_min(j)+(x_max(j)-x_min(j))*rand(1);
   end
   chromo(i,x_num+1:x_num+f_num) = object_fun(chromo(i,:),f_num,x_num,fun);
end

主程式matlab程式碼如下：

%---------------------------------------------------------------------
%程式功能：實現nsga2演算法，測試函式為ZDT1,ZDT2,ZDT3,ZDT4,ZDT6,DTLZ1,DTLZ2
%說明：遺傳運算元為二進位制競賽選擇，模擬二進位制交叉和多項式變異
%      需要設定的引數：pop,gen,pc,pm,yita1,yita2
%作者：(曉風)
%最初建立時間：2018.9.3
%最近修改時間：2018.9.20
%----------------------------------------------------------
clear all
clc
tic;
%引數設定
fun='DTLZ1';
funfun;%函式選擇
pop=300;%種群大小100
gen=250;%250進化代數
pc=1;%交叉概率
pm=1/x_num;%變異概率
% yita1=1;%模擬二進位制交叉引數
% yita2=10;%多項式變異引數
yita1=2;%模擬二進位制交叉引數10
yita2=5;%多項式變異引數50

%%初始化種群
chromo=initialize(pop,f_num,x_num,x_min,x_max,fun);
%%初始種群的非支配排序
[F1,chromo_non]=non_domination_sort(pop,chromo,f_num,x_num);%F為pareto等級為pareto_rank的集合，%p為每個個體p的集合(包括每個個體p的被支配個數n和該個體支配的解的集合s),chromo最後一列加入個體的等級
%%計算擁擠度進行排序
chromo=crowding_distance_sort(F1,chromo_non,f_num,x_num);
i=1;
%%迴圈開始
for i=1:gen
   %%二進位制競賽選擇(k取了pop/2，所以選兩次)
   chromo_parent_1 = tournament_selection(chromo);
   chromo_parent_2 = tournament_selection(chromo);
   chromo_parent=[chromo_parent_1;chromo_parent_2];
   %%模擬二進位制交叉與多項式變異
   chromo_offspring=cross_mutation(chromo_parent,f_num,x_num,x_min,x_max,pc,pm,yita1,yita2,fun );
   %%精英保留策略
   %將父代和子代合併
   [pop_parent,~]=size(chromo);
   [pop_offspring,~]=size(chromo_offspring);
   combine_chromo(1:pop_parent,1:(f_num+x_num))=chromo(:,1:(f_num+x_num));
   combine_chromo((pop_parent+1):(pop_parent+pop_offspring),1:(f_num+x_num))=chromo_offspring(:,1:(f_num+x_num));
   %快速非支配排序
   [pop2,~]=size(combine_chromo);
   [F2,combine_chromo1]=non_domination_sort(pop2,combine_chromo,f_num,x_num);
   %計算擁擠度進行排序
   combine_chromo2=crowding_distance_sort(F2,combine_chromo1,f_num,x_num);
   %精英保留產生下一代種群
   chromo=elitism(pop,combine_chromo2,f_num,x_num);
   if mod(i,10) == 0
       fprintf('%d gen has completed!\n',i);
   end
end
toc;
aaa=toc;

hold on
if(f_num==2)
   plot(chromo(:,x_num+1),chromo(:,x_num+2),'r*');
end
if(f_num==3)
   plot3(chromo(:,x_num+1),chromo(:,x_num+2),chromo(:,x_num+2),'r*');
end
%%--------------------delta度量--------------------------------
% [~,index_f1]=sort(chromo(:,x_num+1));
% d=zeros(length(chromo)-1,1);
% for i=1:(length(chromo)-1)
%     d(i)=((chromo(index_f1(i),x_num+1)-chromo(index_f1(i+1),x_num+1))^2+(chromo(index_f1(i),x_num+2)-chromo(index_f1(i+1),x_num+2))^2)^0.5;
% end
% d_mean1=mean(d);
% d_mean=d_mean1*ones(length(chromo)-1,1);
% d_sum=sum(abs(d-d_mean));
% delta=(d(1)+d(length(chromo)-1)+d_sum)/(d(1)+d(length(chromo)-1)+(length(chromo)-1)*d_mean1);
% disp(delta);
%mean(a)
%(std(a))^2
%--------------------Coverage(C-metric)---------------------
A=PP;B=chromo(:,(x_num+1):(x_num+f_num));%%%%%%%%%%%%%%%%%%%
[temp_A,~]=size(A);
[temp_B,~]=size(B);
number=0;
for i=1:temp_B
   nn=0;
   for j=1:temp_A
       less=0;%當前個體的目標函式值小於多少個體的數目
       equal=0;%當前個體的目標函式值等於多少個體的數目
       for k=1:f_num
           if(B(i,k)<A(j,k))
               less=less+1;
           elseif(B(i,k)==A(j,k))
               equal=equal+1;
           end
       end
       if(less==0 && equal~=f_num)
           nn=nn+1;%被支配個體數目n+1
       end
   end
   if(nn~=0)
       number=number+1;
   end
end
C_AB=number/temp_B;
disp("C_AB:");
disp(C_AB);
%-----Distance from Representatives in the PF(D-metric)-----
A=chromo(:,(x_num+1):(x_num+f_num));P=PP;%%%%%%與論文公式保持相同，注意A和上述不一樣
[temp_A,~]=size(A);
[temp_P,~]=size(P);
min_d=0;
for v=1:temp_P
   d_va=(A-repmat(P(v,:),temp_A,1)).^2;
   min_d=min_d+min(sqrt(sum(d_va,2)));
end
D_AP=(min_d/temp_P);
disp("D_AP:");
disp(D_AP);
filepath=pwd;          
cd('E:\GA\MOEA D\MOEAD_王超\nsga2對比運算元修改\DTLZ1');
save C_AB4.txt C_AB -ASCII
save D_AP4.txt D_AP -ASCII
save solution4.txt chromo -ASCII
save toc4.txt aaa -ASCII
cd(filepath);

3. 實驗結果

本文使用實數編碼的模擬二進位制交叉和多項式變異，交叉概率pc=0.9p_c=0.9，變異概率pm=1/np_m=1/n，ηc=20η_c=20，ηm=20η_m=20。以下為選取的5個非凸非均勻的多目標函式的執行結果如圖4到圖8所示。
函式範圍的matlab程式碼如下：

%測試函式
%--------------------ZDT1--------------------
if strcmp(fun,'ZDT1')
  f_num=2;%目標函式個數
  x_num=30;%決策變數個數
  x_min=zeros(1,x_num);%決策變數的最小值
  x_max=ones(1,x_num);%決策變數的最大值
  load('ZDT1.txt');%匯入正確的前端解
  plot(ZDT1(:,1),ZDT1(:,2),'g*');
  PP=ZDT1;
end
%--------------------ZDT2--------------------
if strcmp(fun,'ZDT2')
  f_num=2;%目標函式個數
  x_num=30;%決策變數個數
  x_min=zeros(1,x_num);%決策變數的最小值
  x_max=ones(1,x_num);%決策變數的最大值
  load('ZDT2.txt');%匯入正確的前端解
  plot(ZDT2(:,1),ZDT2(:,2),'g*');
  PP=ZDT2;
end
%--------------------ZDT3--------------------
if strcmp(fun,'ZDT3')
  f_num=2;%目標函式個數
  x_num=30;%決策變數個數
  x_min=zeros(1,x_num);%決策變數的最小值
  x_max=ones(1,x_num);%決策變數的最大值
  load('ZDT3.txt');%匯入正確的前端解
  plot(ZDT3(:,1),ZDT3(:,2),'g*');
  PP=ZDT3;
end
%--------------------ZDT4--------------------
if strcmp(fun,'ZDT4')
  f_num=2;%目標函式個數
  x_num=10;%決策變數個數
  x_min=[0,-5,-5,-5,-5,-5,-5,-5,-5,-5];%決策變數的最小值
  x_max=[1,5,5,5,5,5,5,5,5,5];%決策變數的最大值
  load('ZDT4.txt');%匯入正確的前端解
  plot(ZDT4(:,1),ZDT4(:,2),'g*');
  PP=ZDT4;
end
%--------------------ZDT6--------------------
if strcmp(fun,'ZDT6')
  f_num=2;%目標函式個數
  x_num=10;%決策變數個數
  x_min=zeros(1,x_num);%決策變數的最小值
  x_max=ones(1,x_num);%決策變數的最大值
  load('ZDT6.txt');%匯入正確的前端解
  plot(ZDT6(:,1),ZDT6(:,2),'g*');
  PP=ZDT6;
end
%--------------------------------------------
%--------------------DTLZ1--------------------
if strcmp(fun,'DTLZ1')
  f_num=3;%目標函式個數
  x_num=10;%決策變數個數
  x_min=zeros(1,x_num);%決策變數的最小值
  x_max=ones(1,x_num);%決策變數的最大值
  load('DTLZ1.txt');%匯入正確的前端解
%     plot3(DTLZ1(:,1),DTLZ1(:,2),DTLZ1(:,3),'g*');
  PP=DTLZ1;
end
%--------------------------------------------
%--------------------DTLZ2--------------------
if strcmp(fun,'DTLZ2')
  f_num=3;%目標函式個數
  x_num=10;%決策變數個數
  x_min=zeros(1,x_num);%決策變數的最小值
  x_max=ones(1,x_num);%決策變數的最大值
  load('DTLZ2.txt');%匯入正確的前端解
%     plot3(DTLZ2(:,1),DTLZ2(:,2),DTLZ2(:,3),'g*');
  PP=DTLZ2;
end
%--------------------------------------------

函式的matlab程式碼如下：

function f = object_fun( x,f_num,x_num,fun )
% --------------------ZDT1--------------------
if strcmp(fun,'ZDT1')
  f=[];
  f(1)=x(1);
  sum=0;
  for i=2:x_num
      sum = sum+x(i);
  end
  g=1+9*(sum/(x_num-1));
  f(2)=g*(1-(f(1)/g)^0.5);
end
% --------------------ZDT2--------------------
if strcmp(fun,'ZDT2')
  f=[];
  f(1)=x(1);
  sum=0;
  for i=2:x_num
      sum = sum+x(i);
  end
  g=1+9*(sum/(x_num-1));
  f(2)=g*(1-(f(1)/g)^2);
end
% --------------------ZDT3--------------------
if strcmp(fun,'ZDT3')
  f=[];
  f(1)=x(1);
  sum=0;
  for i=2:x_num
      sum = sum+x(i);
  end
  g=1+9*(sum/(x_num-1));
  f(2)=g*(1-(f(1)/g)^0.5-(f(1)/g)*sin(10*pi*f(1)));
end
% --------------------ZDT4--------------------
if strcmp(fun,'ZDT4')
  f=[];
  f(1)=x(1);
  sum=0;
  for i=2:x_num
      sum = sum+(x(i)^2-10*cos(4*pi*x(i)));
  end
  g=1+9*10+sum;
  f(2)=g*(1-(f(1)/g)^0.5);
end
% --------------------ZDT6--------------------
if strcmp(fun,'ZDT6')
  f=[];
  f(1)=1-(exp(-4*x(1)))*((sin(6*pi*x(1)))^6);
  sum=0;
  for i=2:x_num
      sum = sum+x(i);
  end
  g=1+9*((sum/(x_num-1))^0.25);
  f(2)=g*(1-(f(1)/g)^2);
end
% --------------------------------------------
% --------------------DTLZ1--------------------
if strcmp(fun,'DTLZ1')
  f=[];
  sum=0;
  for i=3:x_num
      sum = sum+((x(i)-0.5)^2-cos(20*pi*(x(i)-0.5)));
  end
  g=100*(x_num-2)+100*sum;
  f(1)=(1+g)*x(1)*x(2);
  f(2)=(1+g)*x(1)*(1-x(2));
  f(3)=(1+g)*(1-x(1));
end
% --------------------------------------------
% --------------------DTLZ2--------------------
if strcmp(fun,'DTLZ2')
  f=[];
  sum=0;
  for i=3:x_num
      sum = sum+(x(i))^2;
  end
  g=sum;
  f(1)=(1+g)*cos(x(1)*pi*0.5)*cos(x(2)*pi*0.5);
  f(2)=(1+g)*cos(x(1)*pi*0.5)*sin(x(2)*pi*0.5);
  f(3)=(1+g)*sin(x(1)*pi*0.5);
end
% --------------------------------------------
end

3.1 ZDT1

f1=x1f_1=x_1
g=1+9((∑i=2nxi)/(n−1))g=1+9((∑_{i=2}^nx_i )/(n-1))
f2=g(1−(f1/g)0.5).n=30,xi∈[0,1]f_2=g(1-(f_1/g)^{0.5} ).n=30,x_i∈[0,1]
選取種群大小為500，迭代次數500次，pc=0.9,pm=1/30,ηc=20,ηm=20,得到結果如圖4(上)所示。
選取種群大小為100，迭代次數2000次，pc=0.9,pm=1/30,ηc=20,ηm=20,得到結果如圖4(下)所示。
在這裡插入圖片描述
在這裡插入圖片描述
圖4 ZDT1 pareto最優解對比圖（綠色為理論值，紅色為實驗值）

3.2 ZDT2

f1=x1f_1=x_1
g=1+9((∑i=2nxi)/(n−1))g=1+9((∑_{i=2}^n x_i )/(n-1))
f2=g(1−(f1/g)2).n=30,xi∈[0,1]f_2=g(1-(f_1/g)^2 ).n=30,x_i∈[0,1]
選取種群大小為500，迭代次數500次，pc=0.9,pm=1/30,ηc=20,ηm=20,得到結果如圖5（上）所示。
選取種群大小為100，迭代次數2000次，pc=0.9,pm=1/30,ηc=20,ηm=20,得到結果如圖5（下）所示。
在這裡插入圖片描述
在這裡插入圖片描述
圖5 ZDT2 pareto最優解對比圖（綠色為理論值，紅色為實驗值）

3.3 ZDT3

f1=x1f_1=x_1
g=1+9((∑i=2nxi)/(n−1))g=1+9((∑_{i=2}^nx_i )/(n-1))
f2=g(1−(f1/g)0.5−(f1/g)sin⁡(10πf1)).n=30,xi∈[0,1]f_2=g(1-(f_1/g)^{0.5}-(f_1/g)sin⁡(10πf_1)).n=30,x_i∈[0,1]
選取種群大小為136，迭代次數500次，pc=0.9,pm=1/30,ηc=20,ηm=20,得到結果如圖6示。
選取種群大小為136，迭代次數2000次，pc=0.9,pm=1/30,ηc=20,ηm=20,得到結果如圖6示。
在這裡插入圖片描述
在這裡插入圖片描述
圖6 ZDT3 pareto最優解對比圖（綠色為理論值，紅色為實驗值）

3.4 ZDT4

f1=x1f_1=x_1
g=1+9∗10+∑i=2n((xi)2−10cos⁡(4πxi))g=1+9*10+∑_{i=2}^n( (x_i )^2-10cos⁡(4πx_i))
f2=g(1−(f1/g)0.5).n=10,x1∈[0,1];x(2…6)∈[−5,5]f_2=g(1-(f_1/g)^{0.5} ).n=10,x_1∈[0,1];x_(2…6)∈[-5,5]
選取種群大小為100，迭代次數500次，pc=0.9,pm=0.1,ηc=20,ηm=10,得到結果如圖7（上）所示。
選取種群大小為100，迭代次數2000次，pc=0.9,pm=0.1,ηc=20,ηm=10,得到結果如圖7（下）所示。
在這裡插入圖片描述
在這裡插入圖片描述
圖7 ZDT4 pareto最優解對比圖（綠色為理論值，紅色為實驗值）

3.5 ZDT6

f1=1−e−4x1sin6(6πx1)f_1=1-e^{-4x_1} sin^6 (6πx_1)
g=1+9((∑i=2nxi)/(n−1))0.25g=1+9((∑_{i=2}^nx_i )/(n-1))^{0.25}
f2=g(1−(f1/g)2).n=10,xi∈[0,1]f_2=g(1-(f_1/g)^2 ).n=10,x_i∈[0,1]
選取種群大小為100，迭代次數500次，pc=0.9,pm=0.1,ηc=20,ηm=20,得到結果如圖8（上）所示。
選取種群大小為100，迭代次數2000次，pc=0.9,pm=0.1,ηc=20,ηm=20,得到結果如圖8（下）所示。
在這裡插入圖片描述
在這裡插入圖片描述
圖8 ZDT6 pareto最優解對比圖（綠色為理論值，紅色為實驗值）
從結果中可以看出，除ZDT4以外，找到的解幾乎全部是pareto前端上的點，並且解在目標空間上分佈十分均勻，該演算法對於非凸非均勻的多目標函式最優解的尋找還是十分有效的。由於ZDT4具有許多區域性pareto最優前沿，為了擺脫這些區域性前沿，需要對決策向量進行大的改變，因此選取ηm=10，但仍離真正的pareto前沿還有很大的差距。

3.6 對比

表1 γ的均值和方差（500代）

專案 ZDT1 ZDT2 ZDT3 ZDT4 ZDT6
均值 0.0532 0.0972 0.1712 20.2304 0.3284
方差 0.0021 0.0051 0.0048 0.6310 0.0109

表2 Δ的均值和方差（500代）

專案 ZDT1 ZDT2 ZDT3 ZDT4 ZDT6
均值 0.4575 0.4615 0.6640 0.5203 0.9538
方差 0.0021 0.0015 0.0030 0.0013 0.0055

表3 γ的均值和方差（2000代）

專案 ZDT1 ZDT2 ZDT3 ZDT4 ZDT6
均值 0.0011 0.0027 0.0225 12.0341 0.2092
方差 0.0002 0.0002 0.0005 0.3320 0.0021

表4 Δ的均值和方差（2000代）

專案 ZDT1 ZDT2 ZDT3 ZDT4 ZDT6
均值 0.4478 0.4337 0.6586 0.5168 0.4529
方差 0.0001 0.0003 0.0006 0.0018 0.0001

根據上述實驗資料，進行五組實驗，計算distance metric γ，diversity metric Δ的均值和方差，迭代500次後得到結果如表格1-2所示。迭代2000次後得到的結果如表格3-4所示。
      表1顯示了NSGA-Ⅱ（實數編碼500代）在5個問題上執行獲得的distance metric γ的均值和方差，實驗結果表明，除了ZDT4和ZDT6以外的問題均取得了良好的收斂效果，除ZDT4以外，其他問題5次執行結果的差異也都很小。
      表2顯示了NSGA-Ⅱ（實數編碼500代）在5個問題上執行獲得的diversity metric Δ的均值和方差，實驗結果表明，NSGA-Ⅱ在上述問題上均具有較強的傳播多樣性，並且5次執行結果的差異都很小。
      像500代的時候一樣保留所有其他引數，但將代數增加至2000代，得到表3的distance metric γ的均值和方差以及表4的diversity metric Δ的均值和方差。將表3表4與表1表2對比可知，5個問題的distance metric γ隨著代數的增加都有改善，其中ZDT3和ZDT4的改善尤為顯著，除此之外，因為在500代時問題的結果已經有了良好的多樣性分佈了，所以代數的增加並沒有對diversity metric Δ產生巨大的影響。
[1]:Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2):182-197.

相關文章