MATLAB神經網絡（3） 遺傳算法優化BP神經網絡——非線性函數擬合

3.1 案例背景

遺傳算法（Genetic Algorithms）是一種模擬自然界遺傳機制和生物進化論而形成的一種並行隨機搜索最優化方法。

其基本要素包括：染色體編碼方法、適應度函數、遺傳操作和運行參數。

非線性函數：$y=x_{1}^{2}+x_{2}^{2}$

3.2 模型建立

3.2.1 算法流程

 

 

遺傳算法優化使用遺傳算法優化BP神經網絡的權值和闊值，種群中的每個個體都包含了一 個網絡所有權值和闊值，個體通過適應度函數計算個體適應度值，遺傳算法通過選擇、交叉和變異操作找到最優適應度值對應個體。神經網絡預測用遺傳算法得到最優個體對網絡初始權值和閾值賦值，網絡經訓練后預測函數輸出。

神經網絡結構：2-5-1

3.2.2 遺傳算法實現

種群初始化

個體編碼方法為實數編碼，每個個體均為一個實數串，由輸入層與隱含層連接權值、隱含層閾值、隱含層與輸出層連接權值以及輸出層閾值4部分組成。個體包含了神經網絡全部權值和闡值，在網絡結構已知的情況下，就可以構成一個確定的神經網絡。

適應度函數

把預測輸出和期望輸出之間的誤差絕對值和$E$作為個體適應度$F$，計算公式為\[F = k\left( {\sum\limits_{i = 1}^n {{\rm{abs}}({y_i} - {o_i})} } \right)\]$k$為系數。

選擇操作

輪盤賭：基於適應度比例的選擇策略，每個個體$i$的選擇概率$p_{i}$為\[\begin{array}{l}
{f_i} = k/{F_i}\\
{p_i} = \frac{{{f_i}}}{{\sum\limits_{j = 1}^N {{f_j}} }}
\end{array}\]式中，$F_{i}$為個體$i$的適應度值，由於適應度值越小越好，所以在個體選擇前對適應度值取倒數，$k$為系數，$N$為種群個體數目。

交叉操作

由於個體采用實數編碼，所以交叉操作方法采用實數交叉法，第$k$個染色體$a_{k}$和第$l$個染色體$a_{l}$在$j$位的交叉操作方法如下：\[\left\{ \begin{array}{l}
{a_{kj}} = {a_{kj}}(1 - b) + {a_{lj}}b\\
{a_{lj}} = {a_{lj}}(1 - b) + {a_{kj}}b
\end{array} \right.\]式中，$b$是[0,1]之間的隨機數。

變異操作

選取第$i$個個體的第$j$個基因進行變異，變異操作方法如下：\[{a_{ij}} = \left\{ \begin{array}{l}
{a_{ij}} + ({a_{ij}} - {a_{\max }}) * f(g)\;\;r > 0.5\\
{a_{ij}} + ({a_{\min }} - {a_{ij}}) * f(g)\;\;r \le 0.5
\end{array} \right.\]式中，$a_{max}$、$a_{min}$為基因$a_{ij}$的上界和下界，$f(g) = {r_2}{(1 - g/{G_{\max }})^2}$，$r_{2}$為一個隨機數，$g$為當前迭代次數，$G_{max}$為最大進化次數，$r$為[0,1]之間的隨機數。

3.3 編程實現

3.3.1 數據加載

%% 基於遺傳算法神經網絡的預測代碼
% 清空環境變量
clc
clear
% 
%% 網絡結構建立
%讀取數據
input=10*randn(2,2000);
output=sum(input.*input);

%訓練數據和預測數據
input_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);

3.3.2 參數預設

%節點個數
inputnum=2;
hiddennum=5;
outputnum=1;

%構建網絡
net=newff(inputn,outputn,hiddennum);

%% 遺傳算法參數初始化
maxgen=20;                         %進化代數，即迭代次數
sizepop=10;                        %種群規模
pcross=[0.2];                       %交叉概率選擇，0和1之間
pmutation=[0.1];                    %變異概率選擇，0和1之間

%節點總數
numsum=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;

lenchrom=ones(1,numsum);        
bound=[-3*ones(numsum,1) 3*ones(numsum,1)];    %數據范圍

3.3.3 種群初始化

%% 種群初始化
individuals=struct('fitness',zeros(1,sizepop), 'chrom',[]);  %將種群信息定義為一個結構體
avgfitness=[];                      %每一代種群的平均適應度
bestfitness=[];                     %每一代種群的最佳適應度
bestchrom=[];                       %適應度最好的染色體
%初始化種群
for i=1:sizepop
    %隨機產生一個種群
    individuals.chrom(i,:)=Code(lenchrom,bound);    %編碼（binary和grey的編碼結果為一個實數，float的編碼結果為一個實數向量）
    x=individuals.chrom(i,:);
    %計算適應度
    individuals.fitness(i)=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn);   %染色體的適應度
end
FitRecord=[];
%找最好的染色體
[bestfitness,bestindex]=min(individuals.fitness);
bestchrom=individuals.chrom(bestindex,:);  %最好的染色體
avgfitness=sum(individuals.fitness)/sizepop; %染色體的平均適應度
% 記錄每一代進化中最好的適應度和平均適應度
trace=[avgfitness bestfitness]; 

編碼函數：在bound邊界范圍內隨機生成個體。

function ret=Code(lenchrom,bound)
%本函數將變量編碼成染色體，用於隨機初始化一個種群
% lenchrom   input : 染色體長度
% bound      input : 變量的取值范圍
% ret        output: 染色體的編碼值
flag=0;
while flag==0
    pick=rand(1,length(lenchrom));
    ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick; %線性插值，編碼結果以實數向量存入ret中
    flag=test(lenchrom,bound,ret);     %檢驗染色體的可行性
end

檢驗函數：必要的時候可以添加檢驗染色體可行性的代碼。

function flag=test(lenchrom,bound,code)
% lenchrom   input : 染色體長度
% bound      input : 變量的取值范圍
% code       output: 染色體的編碼值
x=code; %先解碼
flag=1;  

適應度函數：以編碼代表的初值進行神經網絡訓練，計算出誤差絕對值和作為適應度。

function error = fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn)
%該函數用來計算適應度值
%x          input     個體
%inputnum   input     輸入層節點數
%outputnum  input     隱含層節點數
%net        input     網絡
%inputn     input     訓練輸入數據
%outputn    input     訓練輸出數據

%error      output    個體適應度值

%提取
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.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
 
%網絡權值賦值
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;

%網絡訓練
net=train(net,inputn,outputn);

an=sim(net,inputn);

error=sum(abs(an-outputn));

3.3.4 進化過程

%% 迭代求解最佳初始閥值和權值
% 進化開始
for i=1:maxgen
    % 選擇
    individuals=Select(individuals,sizepop); 
    avgfitness=sum(individuals.fitness)/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);
    [worestfitness,worestindex]=max(individuals.fitness);
    % 代替上一次進化中最好的染色體
    if bestfitness>newbestfitness
        bestfitness=newbestfitness;
        bestchrom=individuals.chrom(newbestindex,:);
    end
%     individuals.chrom(worestindex,:)=bestchrom;
%     individuals.fitness(worestindex)=bestfitness;
    
    avgfitness=sum(individuals.fitness)/sizepop;
    
    trace=[trace;avgfitness bestfitness]; %記錄每一代進化中最好的適應度和平均適應度
    FitRecord=[FitRecord;individuals.fitness];
end

選擇函數：sumf將[0,1]區間划分為sizepop個區間，生成隨機數落在哪個區間就選取對應的個體。

function ret=Select(individuals,sizepop)
% 本函數對每一代種群中的染色體進行選擇，以進行后面的交叉和變異
% individuals input  : 種群信息
% sizepop     input  : 種群規模
% ret         output : 經過選擇后的種群

%根據個體適應度值進行排序
fitness1=10./individuals.fitness;

sumfitness=sum(fitness1);
sumf=fitness1./sumfitness;
index=[]; 
for i=1:sizepop   %轉sizepop次輪盤
    pick=rand;
    while pick==0    
        pick=rand;        
    end
    for j=1:sizepop    
        pick=pick-sumf(j);        
        if pick<0        
            index=[index j];            
            break;  %尋找落入的區間，此次轉輪盤選中了染色體i，注意：在轉sizepop次輪盤的過程中，有可能會重復選擇某些染色體
        end
    end
end
individuals.chrom=individuals.chrom(index,:);
individuals.fitness=individuals.fitness(index);
ret=individuals;

交叉函數：原為v1=chrom(index(1),pos);，改為v1=chrom(index(1),pos:end);。

function ret=Cross(pcross,lenchrom,chrom,sizepop,bound)
%本函數完成交叉操作
% pcorss                input  : 交叉概率
% lenchrom              input  : 染色體的長度
% chrom     input  : 染色體群
% sizepop               input  : 種群規模
% ret                   output : 交叉后的染色體
 for i=1:sizepop  %每一輪for循環中，可能會進行一次交叉操作，染色體是隨機選擇的，交叉位置也是隨機選擇的，%但該輪for循環中是否進行交叉操作則由交叉概率決定（continue控制）
     % 隨機選擇兩個染色體進行交叉
     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)); %隨機選擇進行交叉的位置，即選擇第幾個變量進行交叉，注意：兩個染色體交叉的位置相同
         pick=rand; %交叉開始
         v1=chrom(index(1),pos:end);
         v2=chrom(index(2),pos:end);
         chrom(index(1),pos:end)=pick*v2+(1-pick)*v1;
         chrom(index(2),pos:end)=pick*v1+(1-pick)*v2; %交叉結束
         flag1=test(lenchrom,bound,chrom(index(1),:));  %檢驗染色體1的可行性
         flag2=test(lenchrom,bound,chrom(index(2),:));  %檢驗染色體2的可行性
         if   flag1*flag2==0
             flag=0;
         else flag=1;
         end    %如果兩個染色體不是都可行，則重新交叉
     end
 end
ret=chrom;

變異函數：

function ret=Mutation(pmutation,lenchrom,chrom,sizepop,num,maxgen,bound)
% 本函數完成變異操作
% pcorss                input  : 變異概率
% lenchrom              input  : 染色體長度
% chrom     input  : 染色體群
% sizepop               input  : 種群規模
% opts                  input  : 變異方法的選擇
% pop                   input  : 當前種群的進化代數和最大的進化代數信息
% bound                 input  : 每個個體的上屆和下屆
% maxgen                input  ：最大迭代次數
% num                   input  : 當前迭代次數
% ret                   output : 變異后的染色體

for i=1:sizepop   %每一輪for循環中，可能會進行一次變異操作，染色體是隨機選擇的，變異位置也是隨機選擇的，
    %但該輪for循環中是否進行變異操作則由變異概率決定（continue控制）
    % 隨機選擇一個染色體進行變異
    pick=rand;
    while pick==0
        pick=rand;
    end
    index=ceil(pick*sizepop);
    % 變異概率決定該輪循環是否進行變異
    pick=rand;
    if pick>pmutation
        continue;
    end
    flag=0;
    while flag==0
        % 變異位置
        pick=rand;
        while pick==0      
            pick=rand;
        end
        pos=ceil(pick*sum(lenchrom));  %隨機選擇了染色體變異的位置，即選擇了第pos個變量進行變異
    
        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,:));     %檢驗染色體的可行性
    end
end
ret=chrom;

3.3.5 結果分析

%% 遺傳算法結果分析 
figure(1)
[r c]=size(trace);
plot([1:r]',trace(:,1),'r');
hold on
plot([1:r]',trace(:,2),'b--');
hold off
title(['適應度曲線  ' '終止代數＝' num2str(maxgen)]);
xlabel('進化代數');ylabel('適應度');
legend('平均適應度','最佳適應度');

 

3.3.6 結果預測

%% 把最優初始閥值權值賦予網絡預測
% %用遺傳算法優化的BP網絡進行值預測
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.trainParam.epochs=100;
net.trainParam.lr=0.1;
%net.trainParam.goal=0.00001;

%網絡訓練
[net,per2]=train(net,inputn,outputn);

%% BP網絡預測
%數據歸一化
inputn_test=mapminmax('apply',input_test,inputps);
an=sim(net,inputn_test);
test_simu=mapminmax('reverse',an,outputps);
error=test_simu-output_test;
plot(error)
title('誤差');
xlabel('測試個體');ylabel('誤差');

精度得到了一定的提高。

3.4 擴展

3.4.1 其他優化方法

粒子群算法、蟻群算法等同樣可以。

3.4.2 網絡結構優化

可以優化隱含層節點數目。

3.4.3 算法的局限性

它只能有限提高原有BP神經網絡的預測精度，並不能把預測誤差較大的BP神經網絡優化為能夠准確預測的BP神經網絡。

其實遺傳算法用處不大，主要還是靠神經網絡，該用神經網絡就用就完了！