[Evolutionary Algorithm] 進化算法簡介

昵稱10504424 2015-12-14

展開全文

閱讀目錄

1. 遺傳算法（Genetic Algorithm，GA）
2. 文化基因算法（Memetic Algorithm，MA）
3. 進化多目標優(yōu)化算法（Multi-Objective Evolutionary Algorithm，MOEA）
4. 參考文獻

　　進化算法，也被成為是演化算法（evolutionary algorithms，簡稱EAs），它不是一個具體的算法，而是一個“算法簇”。進化算法的產生的靈感借鑒了大自然中生物的進化操作，它一般包括基因編碼，種群初始化，交叉變異算子，經營保留機制等基本操作。與傳統(tǒng)的基于微積分的方法和窮舉方法等優(yōu)化算法（具體介紹見博客[Math] 常見的幾種最優(yōu)化方法中的其他數學優(yōu)化方法）相比，進化計算是一種成熟的具有高魯棒性和廣泛適用性的全局優(yōu)化方法，具有自組織、自適應、自學習的特性，能夠不受問題性質的限制，有效地處理傳統(tǒng)優(yōu)化算法難以解決的復雜問題（比如NP難優(yōu)化問題）。

　　除了上述優(yōu)點以外，進化算法還經常被用到多目標問題的優(yōu)化求解中來，我們一般稱這類進化算法為進化多目標優(yōu)化算法（MOEAs）。目前進化計算的相關算法已經被廣泛用于參數優(yōu)化、工業(yè)調度、資源分配、復雜網絡分析等領域。

回到頂部

1. 遺傳算法（Genetic Algorithm，GA）

　　遺傳算法（Genetic Algorithm，簡稱GA）是一種最基本的進化算法，它是模擬達爾文生物進化理論的一種優(yōu)化模型，最早由J.Holland教授于1975年提出。遺傳算法中種群分每個個體都是解空間上的一個可行解，通過模擬生物的進化過程，從而在解空間內搜索最優(yōu)解。

　　遺傳算法的基本操作可以用下圖來描述：

　　個體的編碼方式確定以后，針對上圖操作的具體描述如下：

　　Step 1 種群初始化：根據問題特性設計合適的初始化操作（初始化操作應盡量簡單，時間復雜度不易過高）對種群中的N個個體進行初始化操作；

　　Step 2 個體評價：根據優(yōu)化的目標函數計算種群中個體的適應值（fitness value）；

　　Step 3 迭代設置：設置種群最大迭代次數g_max，并令當前迭代次數g=1；

　　Step 4 個體選擇：設計合適的選擇算子來對種群P(g)個體進行選擇，被選擇的個體將進入交配池中組成父代種群FP(g)，用于交叉變換以產生新的個體。選擇策略要基于個體適應值來進行，假如要優(yōu)化的問題為最小化問題，那么具有較小適應值的個體被選擇的概率相應應該大一些。常用的選擇策略有輪盤賭選擇，錦標賽選擇等。

　　Step 5 交叉算子：根據交叉概率p_m（預先指定，一般為0.9）來判斷父代個體是否需要進行交叉操作。交叉算子要根據被優(yōu)化問題的特性來設計，它是整個遺傳算法的核心，它被設計的好壞將直接決定整個算法性能的優(yōu)劣。

　　Step 6 變異算子：根據變異概率p_c（預先指定，一般為0.1）來判斷父代個體是否需要進行變異操作。變異算子的主要作用是保持種群的多樣性，防止種群陷入局部最優(yōu)，所以其一般被設計為一種隨機變換。

　　通過交叉變異操作以后父代種群FP(g)生成了新的子代種群P(g+1)，令種群迭代次數g=g+1，進行下一輪的迭代操作（跳轉到Step 4），直至迭代次數達到最大的迭代次數。

　　為了更形象說明交叉操作的作用，我們以下圖為例來深入理解一下交叉操作的作用：

　　通過交叉操作，原始的兩個個體組合生成了兩個新的個體組合，這就相當于在解空間進行搜索，每個個體都是解空間的一個可行解。

　　遺傳算法matlab代碼如下：

　　主函數main.m

% This example is used to explain the GA.
%  max f(x1, x2) = 21.5 + x1·sin(4pi*x1) + x2·sin(20pi*x2)
%  s.t.    -3.0 <=  x1 <= 12.1
%           4.1 <=  x2 <= 5.8
clc; clear all;
M = 40;N = 33;
generation = 1000;
it = 1;
pm = 0.05;%變異概率
pc = 0.8;%交叉概率
maxdata = -200;
pop = round(rand(M,N));%初始化矩陣，產生初始種群                      
x1 = decode_x1(pop(:,1:18));  %對x1 x2進行解碼
x2 = decode_x2(pop(:,19:end));
fitness = 21.5 + x1.*sin(4*pi*x1) +x2.*sin(20*pi*x2);%適應度函數
while it < generation
    [B] = seclection(fitness,pop);%輪盤賭選擇
    [newpop] = crossover(pc,B);%交叉
    [B] = mutation(pm,newpop);%變異操作
    pop = B;
    x1 = decode_x1(pop(:,1:18));
    x2 = decode_x2(pop(:,19:end));
    fitness = 21.5 + x1.*sin(4*pi*x1) +x2.*sin(20*pi*x2);%計算適應度
    [fit_best,index] = max(fitness);%本代中的最優(yōu)目標值
    if fit_best >= maxdata
        maxdata = fit_best;
        verter = pop(index,:);
        x_1 = decode_x1(verter(1:18));
        x_2 = decode_x2(verter(19:end));
    end
    num(it) = maxdata;
    it = it + 1;
end
disp(sprintf('max_f=：%.6f',num(it-1)));%輸出最優(yōu)解
disp(sprintf('x_1=：%.5f  x_2=：%.5f',x_1,x_2));%最優(yōu)解對應的x1 x2的值
figure(1)   
plot(num,'k');%繪制圖形

　　變量x1的編碼函數：decode_x1.m

function x1= encode(code)
[M,N] = size(code);
sum = zeros(N);
for i=1:M
    for j=1:N
        sum(i)=sum(i)+code(i,j)*2^(N-j);
    end
    x1(i) = -3.0 + sum(i)*(12.1-(-3.0))/(2^N - 1);
end

　　變量x2的編碼函數：decode_x2.m

function x2=encode_x2(code)
[M,N] = size(code);
sum = zeros(N);
for i=1:M
    for j=1:N
        sum(i)=sum(i)+code(i,j)*2^(N-j);
    end
    x2(i) = 4.1 + sum(i)*(5.8 - 4.1)/(2^N - 1);
end

　　輪盤賭選擇函數：seclection.m

%輪盤賭選擇
function [B]=seclection(fitness,A)
    [M,N]=size(A);
    fit_sum = sum(fitness);
    for i = 1:M
        probability(i) = fitness(i)/fit_sum;
    end
    for i = 2:M
        probability(i) =probability(i) + probability(i-1);
    end
    for j = 1:M
        p = rand();
        i = 1;
        while p > probability(i)
              i = i+1;
        end
        B(j,:) = A(i,:);
    end

　　個體交叉算子：crossover.m

%交叉
function [newpop]=crossover(pc,A)
newpop=A;
[M,N]=size(A);
    for i= 1:2:M-1
        if rand < pc
            p = ceil(rand()*N);
            for j= p:N
                ch = A(i,j);
                newpop(i,j) = A(i+1,j);
                newpop(i+1,j) = ch;
            end
        end
    end

　　個體變異算子：mutation.m

%變異
function [newpop]=mutation(pm,A)
[M,N]=size(A);
newpop=A;
W = rand(M,N)<pm;
    for i=1:M
        for j=1:N
            if W(i,j) ==1
               if A(i,j)~=1
                  newpop(i,j)= 1;
               else
                  newpop(i,j) = 0;
               end
           end
       end
    end

回到頂部

2. 文化基因算法（Memetic Algorithm，MA）

　　文化基因算法（Memetic Algorithm，簡稱MA），也被稱為是“密母算法”，它是由Mpscato在1989年提出的。文化基因算法是一種基于種群的全局搜索和基于個體的局部啟發(fā)式搜索的結合體，它的本質可以理解為：

Memetic = GA + Local Search

即memetic算法實質上為遺傳算法加上一個局部搜索（Local Search）算子。局部搜索算子可以根據不同的策略進行設計，比如常用的爬山機制、模擬退火、貪婪機制、禁忌搜索等。

回到頂部

3. 進化多目標優(yōu)化算法（Multi-Objective Evolutionary Algorithm，MOEA）

　　對于一個優(yōu)化問題而言，如果其只有一個目標方程，那么我們稱之為單目標優(yōu)化問題；而一旦方程個數達到兩個或者兩個以上，那么它被相應的稱之為多目標優(yōu)化問題（Multi-objective Optimization Problems，簡稱為MOPs）。

　　對于一個多目標優(yōu)化問題而言，問題的最優(yōu)解可能不止一個，而應該是一組，我們通常稱這組最優(yōu)解為相應多目標優(yōu)化問題的一個非支配解集，或者稱為是Pareto解集，其中的每一個解我們稱之為Pareto解（Pareto是引自一個經濟學的術語）。求解多目標優(yōu)化問題的解法有很多，比如常見的目標規(guī)劃方法，目標分解方法，目標化多為少方法（將多個目標表示為一個）等。進化算法在解決多目標問題上有著天然的優(yōu)勢，對于一個進化多目標優(yōu)化算法而言，它可以對多個目標函數同時進行優(yōu)化，而且輸出一組非支配的Pareto解集，從而可以有效地求解多目標問題。

　　下圖為一個具有兩目標的多目標優(yōu)化問題的Pareto解集示意圖：

　　常用的進化多目標優(yōu)化算法有Deb提出基于擁擠度距離度量的NSGA-II，張青富老師提出的基于分解思想的MOEA/D算法，以及西電公老師提出的基于免疫克隆機制的NNIA算法等。最近，Deb等人在NSGA-II的基礎上又提出一種針對2~16個目標函數同時優(yōu)化的高維問題的NSGA-III算法，對進化多目標優(yōu)化感興趣的博友可以深入的看看相關的參考文獻，這方面的內容還是很有意思的。

　　NSGA-II算法的matlab代碼參考如下：

　　主函數main.m

clear all
clc
tic
N=50;%the number of the population
numvariate=1;%the number of the variate of each individual
numfun=2;
minx=-10^3;
maxx=10^3;
iteration=1;
Pt0=initializationPt0(minx,maxx,N,numvariate,numfun);

while(iteration<500)
fprintf('-->iteration= %d\n',iteration);
[Rt]=Genetic_Operators( Pt0,minx,maxx,numvariate,numfun);%%執(zhí)行交叉變異后，子代個體的支配面號，擁擠度距離，這兩項還沒計算
[ Fi] = fast_nondominated_sort( Rt,numvariate,numfun);%Fi存儲占優(yōu)面上相應每個個體在種群Rt中的編號，即行號
for i=1:2*N
    numFi=Fi(i,1);
    if (numFi~=0)
         [Idistance]=crowdiing_distance_assignment(Fi(i,:),Rt,numvariate,numfun);
         for j=1:numFi
            Rt(Fi(i,j+1),numvariate+numfun+1)=i;
             Rt(Fi(i,j+1),numvariate+numfun+2)=Idistance(j);
         end
    end
end
Rt;

Pt1=zeros(1,N+1);
i=1;
 while((Pt1(1)+Fi(i,1))<=N)
   Idistance1=crowdiing_distance_assignment(Fi(i,:),Rt,numvariate,numfun);
   for j=1:Fi(i,1)
        Pt1(1)=Pt1(1)+1;
       Pt1(Pt1(1)+1)=Fi(i,j+1);
   end
   i=i+1;
 end
if Pt1(1)<N
    Idistance=crowdiing_distance_assignment(Fi(i,:),Rt,numvariate,numfun);
   [valueIdis,rankIdis]=sort(Idistance);
   NumFi=Fi(i,1);
    for j=(Pt1(1)+2):(N+1)
        k=j-Pt1(1);
        Pt1(j)=Fi(i,rankIdis(NumFi-k+2)+1);
    end
end
b=Pt1(2:N+1);
Pt0=Rt(b,:);
iteration=iteration+1;
end
f1=Pt0(:,(numvariate+1));
f2=Pt0(:,(numvariate+2));
figure;
plot(f1,f2,'bo');
xlabel('f1');
ylabel('f2');
toc

　　目標函數：funfvalue.m

%%%input :the population Rt of size N*m,output:the value of the population
%%%Rt of size N*numf,numf if the number of the objective 
function [ funvalueRt ] = funfvalue( Rt)
%% 測試數據SCH
  funvalueRt(:,1)=Rt.^2;
  funvalueRt(:,2)=(Rt-2).^2;

　　種群初始化函數：GeneratePt0.m

% %The interval [a, b] is divided into N/10 equal parts,
function [ Pt0 ] = GeneratePt0(minx,maxx,N,numvariate,numfun)
Pt0=unifrnd(minx,maxx,N,numvariate);%種群數量為100
fPt0=zeros(N,numfun+2);
Pt0=cat(2,Pt0,fPt0);
end

　　基因操作（交叉和變異）：Genetic_Operators.m

% % input Pt ,size of N*m;output Rt,size of 2N*m
function [ Rt ] = Genetic_Operators( Pt,minx,maxx,numvariate,numfun)
%% 不交叉則變異，避免種群中數值相同的個體存在，這樣可以保持種群的多樣性
pc=0.9;
pm=0.1;
cetac=20;
cetam=20;
maxminx=maxx-minx;
N=size(Pt,1);
numvariate;
Qt=Pt;
crossPt=zeros(N,numvariate);
for i=1:numvariate
crossPt(:,i)=randperm(N);
end
detamax=10;%這個參數有待調整
for j=1:numvariate %%先交叉
    for i=1:ceil(N/2)
        pccross=rand(1);
         if(pccross<pc)
             ui=rand(1);
             if ui<=0.5
                 betau=(2*ui)^(1/(cetac+1));
             else 
                 if ui==1
                     ui=1-0.005;
                 end
                 betau=(1/(2*(1-ui)))^(1/(cetac+1));
             end
            Qt1=0.5*((1-betau)*Pt(crossPt(i,j),j)+(1+betau)*Pt(N-crossPt(i,j)+1,j));
            Qt2=0.5*((1+betau)*Pt(crossPt(i,j),j)+(1-betau)*Pt(N-crossPt(i,j)+1,j));
            Qt(2*i-1,j)=mod(Qt1-minx,maxminx)+minx;
            Qt(2*i,j)=mod(Qt2-minx,maxminx)+minx;
         else
             pmu1=rand(1);
             if pmu1<0.5
                deta1=(2*pmu1)^(1/(cetam+1))-1;
             else deta1=1-(2*(1-pmu1))^(cetam+1);
             end
             pmu2=rand(1);
             if pmu2<0.5
                deta2=(2*pmu2)^(1/(cetam+1))-1;
             else deta2=1-(2*(1-pmu2))^(cetam+1);
             end
            Qt1=Pt(crossPt(i,j),j)+deta1*detamax;
            Qt2=Pt(N-crossPt(i,j)+1,j)+deta2*detamax;
            Qt(2*i-1,j)=mod(Qt1-minx,maxminx)+minx;
            Qt(2*i,j)=mod(Qt2-minx,maxminx)+minx;
          end
     end
end

funvalueQt=funfvalue(Qt(:,1:numvariate));
Qt(:,(numvariate+1):(numvariate+numfun))=funvalueQt;
Rt=cat(1,Pt,Qt);
end

　　快速的非支配排序函數：fast_nondominated_sort.m

% %%Fii存儲占優(yōu)面上相應每個個體在種群Rt中的編號，即行號
function [Fii] = fast_nondominated_sort( Rt,numvariate,numfun)
N=size(Rt,1);
Sp=zeros(N,N+1);%Sp的第一列表示每個Spi所擁有的個數，即p所支配的個數，從第二列開始表示被p支配的個體
np=zeros(N,1);
Finum=zeros(N,N+1);
% Fi=zeros(N,N+1,numvariate);
prank=zeros(N,1);
funvalueRt=Rt(:,(numvariate+1):(numvariate+numfun));
for p=1:N
    for q=1:N
        if q~=p
            %%% 這里還有另外一種情況，即Rt(p)不支配Rt(q)，也不被Rt(q)支配
            if dominate(funvalueRt(p,:),funvalueRt(q,:))==1
                Sp(p,1)=Sp(p,1)+1;
                Sp(p,Sp(p,1)+1)=q;
            elseif dominate(funvalueRt(q,:),funvalueRt(p,:))==1
                    np(p)=np(p)+1;
            end
        end
    end
    
    if np(p)==0
        prank(p)=1;
        Finum(1,1)=Finum(1,1)+1;
        Finum(1,Finum(1,1)+1)=p;
    end
end
Sp;
% np;
i=1;
while Finum(i,1)~=0
    Q=[];
%     fprintf('-->i= %d\n',i);
    Finum(i,:);
    for ip=2:(Finum(i,1)+1)
        p=Finum(i,ip);
        for iq=2:(Sp(p,1)+1)
            q=Sp(p,iq);
            np(q)=np(q)-1;
            if np(q)==0
                prank(q)=i+1;
                Q=[Q,q];
            end
        end
    end
    i=i+1;
    nun2Q=size(Q,2);
    Finum(i,1)=nun2Q;
    for j=1:nun2Q
        Finum(i,j+1)=Q(j);
    end
end
Fii=Finum;
end

　　擁擠度距離計算函數：

%%% input:Fi of size 1*(2*N+1),output:Idistance of size 1*(FiRt(1,1))
function [Idistance] = crowdiing_distance_assignment(Fi,Rt,numvariate,numfun)
NFi=Fi(1,1);
Idistance=zeros(1,NFi);
funvalueI1=zeros(NFi,numfun);
if NFi<=2
    Idistance=Inf*ones(1,NFi);
else 
    b=Fi(2:(Fi(1,1)+1));
    funvalueI1=Rt(b,(numvariate+1):(numvariate+numfun));
%     b=Fi(2:(Fi(1,1)+1));
%     for i=1:NFi
%         funvalueI1(i,:)=Rt(b(i),(numvariate+1):(numvariate+numfun));
%     end
    
%   fprintf('-->Idistancel= %d\n',l);
    for i=1:numfun
        [valueI1,rankI1]=sort(funvalueI1(:,i));
        Idistance(rankI1(1))=Inf;
        Idistance(rankI1(NFi))=Inf;
        for j=2:(NFi-1)
            Idistance(rankI1(j))=Idistance(rankI1(j))+(valueI1(j+1)-valueI1(j-1))/(valueI1(NFi)-valueI1(1));
        end
    end
end
end

　　支配函數：dominate.m

%% 如果Rtp支配Rtq，返回值value=1，否則value=0；
function [value] = dominate(Rtp,Rtq)
funvaluepq=Rtp-Rtq;
value=0;
if funvaluepq<0
    value=1;
elseif funvaluepq<=0
    if(sum(funvaluepq)<0)
        value=1;
    else value=0;
    end
end
end