進化算法,也被成為是演化算法(evolutionary algorithms,簡稱EAs),它不是一個具體的算法,而是一個“算法簇”。進化算法的產生的靈感借鑒了大自然中生物的進化操作,它一般包括基因編碼,種群初始化,交叉變異算子,經營保留機制等基本操作。與傳統的基於微積分的方法和窮舉方法等優化算法(具體介紹見博客[Math] 常見的幾種最優化方法中的其他數學優化方法)相比,進化計算是一種成熟的具有高魯棒性和廣泛適用性的全局優化方法,具有自組織、自適應、自學習的特性,能夠不受問題性質的限制,有效地處理傳統優化算法難以解決的復雜問題(比如NP難優化問題)。
除了上述優點以外,進化算法還經常被用到多目標問題的優化求解中來,我們一般稱這類進化算法為進化多目標優化算法(MOEAs)。目前進化計算的相關算法已經被廣泛用於參數優化、工業調度、資源分配、復雜網絡分析等領域。
1. 遺傳算法(Genetic Algorithm,GA)
遺傳算法(Genetic Algorithm,簡稱GA)是一種最基本的進化算法,它是模擬達爾文生物進化理論的一種優化模型,最早由J.Holland教授於1975年提出。遺傳算法中種群分每個個體都是解空間上的一個可行解,通過模擬生物的進化過程,從而在解空間內搜索最優解。
遺傳算法的基本操作可以用下圖來描述:
個體的編碼方式確定以后,針對上圖操作的具體描述如下:
Step 1 種群初始化:根據問題特性設計合適的初始化操作(初始化操作應盡量簡單,時間復雜度不易過高)對種群中的N個個體進行初始化操作;
Step 2 個體評價:根據優化的目標函數計算種群中個體的適應值(fitness value);
Step 3 迭代設置:設置種群最大迭代次數gmax,並令當前迭代次數g=1;
Step 4 個體選擇:設計合適的選擇算子來對種群P(g)個體進行選擇,被選擇的個體將進入交配池中組成父代種群FP(g),用於交叉變換以產生新的個體。選擇策略要基於個體適應值來進行,假如要優化的問題為最小化問題,那么具有較小適應值的個體被選擇的概率相應應該大一些。常用的選擇策略有輪盤賭選擇,錦標賽選擇等。
Step 5 交叉算子:根據交叉概率pm(預先指定,一般為0.9)來判斷父代個體是否需要進行交叉操作。交叉算子要根據被優化問題的特性來設計,它是整個遺傳算法的核心,它被設計的好壞將直接決定整個算法性能的優劣。
Step 6 變異算子:根據變異概率pc(預先指定,一般為0.1)來判斷父代個體是否需要進行變異操作。變異算子的主要作用是保持種群的多樣性,防止種群陷入局部最優,所以其一般被設計為一種隨機變換。
通過交叉變異操作以后父代種群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);%本代中的最優目標值 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)));%輸出最優解 disp(sprintf('x_1=:%.5f x_2=:%.5f',x_1,x_2));%最優解對應的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年提出的。文化基因算法是一種基於種群的全局搜索和基於個體的局部啟發式搜索的結合體,它的本質可以理解為:
Memetic = GA + Local Search
即memetic算法實質上為遺傳算法加上一個局部搜索(Local Search)算子。局部搜索算子可以根據不同的策略進行設計,比如常用的爬山機制、模擬退火、貪婪機制、禁忌搜索等。
3. 進化多目標優化算法(Multi-Objective Evolutionary Algorithm,MOEA)
對於一個優化問題而言,如果其只有一個目標方程,那么我們稱之為單目標優化問題;而一旦方程個數達到兩個或者兩個以上,那么它被相應的稱之為多目標優化問題(Multi-objective Optimization Problems,簡稱為MOPs)。
對於一個多目標優化問題而言,問題的最優解可能不止一個,而應該是一組,我們通常稱這組最優解為相應多目標優化問題的一個非支配解集,或者稱為是Pareto解集,其中的每一個解我們稱之為Pareto解(Pareto是引自一個經濟學的術語)。求解多目標優化問題的解法有很多,比如常見的目標規划方法,目標分解方法,目標化多為少方法(將多個目標表示為一個)等。進化算法在解決多目標問題上有着天然的優勢,對於一個進化多目標優化算法而言,它可以對多個目標函數同時進行優化,而且輸出一組非支配的Pareto解集,從而可以有效地求解多目標問題。
下圖為一個具有兩目標的多目標優化問題的Pareto解集示意圖:
常用的進化多目標優化算法有Deb提出基於擁擠度距離度量的NSGA-II,張青富老師提出的基於分解思想的MOEA/D算法,以及西電公老師提出的基於免疫克隆機制的NNIA算法等。最近,Deb等人在NSGA-II的基礎上又提出一種針對2~16個目標函數同時優化的高維問題的NSGA-III算法,對進化多目標優化感興趣的博友可以深入的看看相關的參考文獻,這方面的內容還是很有意思的。
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);%%執行交叉變異后,子代個體的支配面號,擁擠度距離,這兩項還沒計算 [ Fi] = fast_nondominated_sort( Rt,numvariate,numfun);%Fi存儲占優面上相應每個個體在種群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存儲占優面上相應每個個體在種群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
4. 參考文獻
[1] 進化多目標優化算法研究,軟件學報,2009.
[2] Messy genetic algirithms: Motivation, analysis, and first results. Complex systems, 1989.
[3] Genetic algorithms in search, optimization, and mechine learning. Addion wesley, 1989.
[4] A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transaction on Evolutionary Computation, 2002.