MATLAB神經網路（三）：遺傳演算法優化BP

01-26

遺傳演算法是1962年由Holland提出的模擬自然界遺傳機制和生物進化論而成的一泓並行隨機搜索最優化方法。它把自然界優勝劣汰、適者生存的生物進化原理引入優化參數形成的編碼串聯群體中，按照所選擇的適應度函數並通過遺傳中的選擇、交叉和變異對個體進行篩選，使適應度值好的個體被保留，適應度較差的個體被淘汰，新的群體既繼承了上一代的信息，又優於上一代。這樣反覆循環，直至滿足條件。遺傳演算法基本的操作為：

1.選擇操作：從舊群體中以一定概率選擇個體到新群體中，個體被選中的概率跟適應度值有關，個體適應度值越好，被選中的概率越大。

2。交叉操作：從個體中選擇兩個個體，通過兩個染色體的交換組合，來產生新的優秀個體。交叉過程為從群體中任選兩個染色體，隨機選擇一點或多點染色體位置進行交換。交叉操作如圖。

3.變異操作：從群體中任選一個個體，選擇染色體中的一點進行變異以產生更優秀的個體。

遺傳演算法具有高效啟發式搜索、並行計算等特點，目前已經應用在函數優化、組合優化以及生產調度等方面。

遺傳演算法的基本要素：染色體編碼方法、適應度函數、遺傳操作和運行參數。染色體編碼方法是指個體的編碼方法，目前包括二進位法、實數法等。二進位法把個體編碼成為一個二進位串，實數法是指把個體編碼成為一個實數串。適應度函數是指根據進化目標編寫的計算個體適應度值的函數。遺傳操作是指選擇操作、交叉操作和變異操作。

運行參數是演算法在初始化時確定的參數，主要包括群體大小M，遺傳代數G，交叉概率P和變異概率Pm。

我們通過遺傳演算法優化最基礎的BP神經網路，以擬合非線性函數：

$y^2=x_1^2+x_2^2$

本例中，由於擬合非線性函數有2個輸入參數，1個輸出參數，所以設置的BP神經網路結構為2-5-1，共有2*5+5*1個權值，5+1=6個閾值，所以遺傳演算法個體編碼長度為16+5=21。

適應度函數用訓練數據訓練BP神經網路，並把訓練數據預測誤差作為個體適應度。

function error=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn)%該函數用來計算適應度值%x input個體%inputnum input 輸入層節點數%hiddennum input 隱含層節點數%outputnum input 輸出層節點數%net input 網路%inputn input 訓練輸入數據%outputn input 訓練輸出數據w1=x(1:inputnum*hiddennum);B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum)B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum)net.iw(1,1)=reshape(w1,hiddennum,inputnum);net.lw(2,1)=reshape(w2,outputnum,hiddennum);net.b(1)=reshape(B1,hiddennum,1);net.b(2)=B2;%BP神經網路構建net=newff(inputn,outputn,hiddennum);net.trainParam.epochs=20;net.trainParam.lr=0.1;net.trainParam.goal=0.00001;net.trainParam.show=100;net.trainParam.showWindow=0;%神經網路訓練net=train(net,inputn,outputn);%BP神經網路預測an=sim(net,inputn);%預測誤差和作為個體適應度值error=sum(abs(an-outputn));

選擇操作採用輪盤賭法從種群中選擇適應度好的個體組成新種群。

function ret=select(individuals,sizepop)%該函數用於進行選擇操作%individuals input 種群信息%sizepop input 種群規模%ret output 選擇後的新種群%求適應度值倒數fitness1=10./individuals.fitness; %individuals.fitness為個體適應度值%個體選擇概率sumfitness=sum(fitness1);sumf=fitness1./sumfitness;%採用輪盤賭法選擇新個體index=[];for i=1:sizepop %sizepop為種群數 pick=rand; while pick==0 pick=rand; end for i=1:sizepop pick=pick-sumf(i); if pick<0 index=[index i]; break; end endend%新種群individuals.chrom=individuals.chrom(index,:); %individuals.chrom為種群中個體individuals.fitness=individuals.fitness(index);ret=individuals;

交叉操作：從種群中選擇兩個個體，按一定概率交叉得到新個體。

function ret=Cross(pcross,lenchrom,chrom,sizepop,bound)%該函數用於進行交叉操作%pcross 交叉概率%lenchrom 個體長度%chrom 種群個體%sizepop 種群規模%ret 交叉後的新種群for i=1:sizepop %隨機選擇兩個個體 pick=rand(1,2); while prod(pick)==0 pick=rand(1,2); end index=ceil(pick,*sizepop)判斷是否交叉 pick=rand; while pick==0 pick=rand; end if pick>pcross continue; end flag=0; while flag==0 %選擇交叉位置 pick=rand; while pick==0 pick=rand; end pos=ceil(pick,*sum(lenchrom)); %lenchrom為個體長度%個體交叉pick=rand;v1=chrom(index(1),pos);v2=chrom(index(2),pos);chrom(index(1),pos)=pick*v2+(1-pick)*v1;chrom(index(2),pos)=pick*v1+(1-pick)*v2;%測試新個體是否滿足約束要求flag1=test(lenchrom,bound,chrom(index(1),:));flag2=test(lenchrom,bound,chrom(index(2),:)); if flag1*flag2==0 flag=0; else flag=1; end end endret=chrom

變異操作：從種群中隨機選擇一個個體，按一定概率變異得到新個體。

function ret=Mutation(pmutation,lenchrom,chrom,sizepop,num,maxgen,bound)%pmutation 變異概率%lenchrom 個體長度%chrom 種群個體%sizepop 種群規模%bound 個體上界和下界%maxgen 最大迭代次數%num 當前迭代次數%ret 交叉後的新種群for i=1:sizepop pick=rand; while pick==0 pick=rand; end index=ceil(pick*sizepop); %判斷是否變異 pick=rand; if pick>pmutation %pmutation 為變異概率 continue; end flag=0; while flag==0 %隨機選擇變異位置 pick=rand; while pick==0 pick=rand; end pos=ceil(pick*sum(lenchrom)); %lenchrom為個體長度 %變異操作 v=chrom(i,pos); v1=v-bound(pos,1); v2=bound(pos,2)-v; pick=rand; fg=(rand*(1-num/maxgen)^2; if pick>0.5 chrom(i,pos)=chrom(i,pos)+(bound(pos,2)-chrom(i,pos))*fg; else chrom(i,pos)=chrom(i,pos)-(chrom(i,pos)-bound(pos,1))*fg; end flag=test(lenchrom,bound,chrom(i,:)); %新個體是否滿足約束要求 endendret=chrom

遺傳演算法主函數

load data input output%讀取數據%網路結構inputnum=2;hiddennum=5;outputnum=1;%取訓練數據和預測數據input_train=input(1:1900,:);input_test=input(1901:2000,:);output_train=output(1:1900);output_test=output(1901:2000);%數據歸一化intput_train=input(1:1900,:);input_test=input(1901:2000,:);output_train=output(1:1900);output_test=output(1901:2000);%數據歸一化[inputn,inputps]=mapminmax(input_train);[outputn,outputps]=mapminmax(output_train);%構建網路net=newff(intputn,outputn,hiddennum);%遺傳演算法參數初始化maxgen=50; %迭代次數sizepop=10; %種群規模pcross=[0.4]; %交叉概率pmutation=[0.2]; %變異概率numsum=inputnum*hiddennum+hiddennum+hiddenum*outputnum+outputnum;lenchrom=cones(1,numsum); %個體長度bound=[-3*ones(numsum,1)3*ones(numsum,1)]; %個體範圍%種群信息定義為結構體indivisuals=struct(fitness,zeros(1,sizepop),chrom,[]);avgfitness=[]; %每代平均適應度值bestfitness=[]; %每代最佳適應度值bestchrom=[]; %最優個體%計算個體適應度值 for i=1:sizepop %個體初始化 individuals.chrom(i,:)=Code(lenchrom,bound); %計算個體適應度值 x=individuals.chrom(i,:); individuals.fitness(i)=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn); end %迭代尋優 for i=1:maxgen %選擇操作 individuals=Select(individuals,sizepop); %交叉操作 individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,bound); %變異操作 individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,i,maxgen,bound); %計算適應度值 for j=1:sizepop x=individuals.chrom(j,:);%個體 individuals.fitness(j)=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn); end %尋找最優最差個體 [newbestfitness,newbestindex]=min(individuals.fitness); [newworstfitness,newworstindex]=max(individuals,fitness); %最優個體更新 if bestfitness>newbestfitness bestfitness=newbestfitness; bestchrom=individuals.chrom(newbestindex,:); end individuals.chrom(newworstindex,:)=bestchrom individuals.fitness(newworstindex)=bestfitness; %記錄最優個體適應度值和平均適應度值avgfitness=sum(individuals.fitness)/sizepop;trace=[trace;avgfitness bestfitness];end%把最優個體賦給BP神經網路，用該網路擬合非線性函數x=bestchromw1=x(1:inputnum*hiddennum);B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);w2=x(inputnum*hiddennum+hiddennum+1:inputnnum*hiddennum+hiddenum+hiddennum*outputnum);B2=x(inputnnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);net.iw(1,1)=reshape(w1,hiddennum,inputnum);net.lw(2,1)=reshape(w2,outputnum,hiddennum);neet.b(1)=reshape(B1,hiddennum,1);net.b(2)=B2;%神經網路參數net.trainParam.epochs=100;net,trainParam.lr=0.1;%net.trainParam.goal=0.00001;%BP神經網路訓練[net,per2]=train(net,inputn,outputn);%BP神經網路預測inputn_test=mapminmax(apply,input)test,inputps);an=sim(net,input_test);%預測結果反歸一化test_simu=mapminmax(reverse,an,outputps);