接我上一條的筆記，今天我加了遺傳算法進行優化BP神經網絡，其基本思想就是用個體代表網絡的初始權值和門檻值、個體值初始化的BP神經網絡的預測誤差作為該個體的适應度值，通過選擇、交叉、變異操作尋找最優個體，即最優的BP神經網絡初始權值。
了解過BP神經網絡的朋友都知道，BP神經網絡初始神經元之間的權值和門檻值一般随機選擇，是以容易陷入局部最小值，通過引入遺傳算法就可以很好找到全局最小值，廢話不多說，直接上。（代碼僅供參考，小白起步）
題目：拟合非線性函數：y=x1x1+x2x2（平方和）
設定網絡：輸入層2個神經元，隐含層5個神經元，輸出層1個神經元
先來了解一下種群、個體、染色體、基因都分别是什麼？種群由全部個體組成，其包含了所有個體的資訊，由此我們可以聯想到可以定義一個結構體來表示種群，來囊括所有個體的資訊；個體，在遺傳算法中我們可以了解成解決一個問題的方案，一個問題有不同種解決方案（也就是不同的個體），模拟自然界中選擇、淘汰個體以得到對解決問題最有用的個體；而一個個體中，有染色體，染色體又由基因組成，是以針對本案例我們不難了解，染色體其實由4部分構成：輸入層到隐含層之間的連接配接權值（一個二維矩陣）、隐含層到輸出層之間的連接配接權值（一個二維矩陣）、隐含層神經元門檻值（一維矩陣）、輸出層神經元門檻值（一維矩陣），而基因就是神經元與神經元之間的連接配接權值（一個數值）。
具體會發生什麼、怎麼發生，我們通過下面代碼一步一步來：
%% 該代碼為基于遺傳算法神經網絡的預測代碼
% 清空環境變量
clc
clear
%
%%網絡結建構立
%讀取資料
load data input output
%節點個數
inputnum=2;
hiddennum=5;
outputnum=1;
%訓練資料和預測資料
input_train=input(1:1900,:)';%取1到1900行的資料來訓練
input_test=input(1901:2000,:)';%取1901到2000行的資料來測試
output_train=output(1:1900);
output_test=output(1901:2000);
%選連樣本輸入輸出資料歸一化
[inputn,inputps]=mapminmax(input_train);%歸一化到[-1,1]之間，inputps用來作下一次同樣的歸一化
[outputn,outputps]=mapminmax(output_train);
%建構網絡
net=newff(inputn,outputn,hiddennum);%單隐含層，5個隐含層神經元
%% 遺傳算法參數初始化
maxgen=20; %進化代數，即疊代次數
sizepop=10; %種群規模
pcross=0.2; %交叉機率選擇，0和1之間
pmutation=0.1; %變異機率選擇，0和1之間
%節點總數：輸入隐含層權值、隐含門檻值、隐含輸出層權值、輸出門檻值（4個基因組成一條染色體）
numsum=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;%21個，10,5,5,1
lenchrom=ones(1,numsum);%個體長度，暫時先了解為染色體長度，是1行numsum列的矩陣
bound=[-3*ones(numsum,1) 3*ones(numsum,1)]; %是numsum行2列的串聯矩陣，第1列是-3，第2列是3
%------------------------------------------------------種群初始化--------------------------------------------------------
individuals=struct('fitness',zeros(1,sizepop), 'chrom',[]); %将種群資訊定義為一個結構體：10個個體的适應度值，10條染色體編碼資訊
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
通過上面代碼我們就知道，結構體individuals就是種群，裡面包含了兩個矩陣，一個用來放全部個體的适應度值，一個用來存放染色體的資訊。好，函數從Code開始，也就是給染色體進行編碼，讓每個基因都指派初始值。
function ret=Code(lenchrom,bound)
%本函數将變量編碼，用于随機初始化一個種群
% lenchrom input : 染色體長度（其實就是1行21列的矩陣）
% bound input : 變量的取值範圍
% ret output: 染色體的編碼值
flag=0;
while flag==0
pick=rand(1,length(lenchrom));%length得到一個矩陣裡較大的行數或者列數，lenchrom是1行numsum列矩陣故傳回numsum，即pick是1行numsum列的随機數矩陣
%bound(:,1)'為取1行numsum列的值都是-3，bound(:,2)為取numsum行第2列的值都是3，(bound(:,2)-bound(:,1))'得到1行numsum列矩陣為6再與pick逐個元素相乘
ret=bound(:,1)'+(bound(:,2)-bound(:,1))'.*pick; %線性插值，編碼結果以實數向量存入ret中，ret是1行numsum列的矩陣
flag=test(ret); %檢驗染色體的可行性
end
編碼之後，每個基因都有了随機數初始值，我們再接着看test函數裡有什麼。
function flag=test(code)
%判斷門檻值和權值是否超界
% lenchrom input : 染色體長度
% bound input : 變量的取值範圍
% code output: 染色體的編碼值
% x=code; %先解碼，x是1行numsum列的矩陣
% flag=1;
f1=isempty(find(code>3, 1));%isempty(A)若A為空，則傳回1，否則傳回0（這個3其實是bound的邊界值）
f2=isempty(find(code<-3, 1));%find中如果滿足條件，則傳回元素位置，最後面的1是為了讓程式提高尋找效率
%fine(code<-3)能夠找出所有滿足條件的元素位置，但我們隻需要檢測code<-3就好了沒必要找到全部，是以用find(code<-3,1)
if f1*f2==0 %有基因（權值或者門檻值）超界，傳回0讓程式重新編碼
flag=0;
else
flag=1;%編碼符合條件，傳回1跳出編碼
end
% if (x(1)<0)&&(x(2)<0)&&(x(3)<0)&&(x(1)>bound(1,2))&&(x(2)>bound(2,2))&&(x(3)>bound(3,2))
% flag=0;
% end
接着就是我們的fun函數（要拉上去再看看啦），fun函數是用來給每個個體生成适應度值的，初始值是随機數
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);%取到輸入層與隐含層連接配接的權值，在Code函數裡已經指派
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);%将w1由1行inputnum*hiddennum列轉為hiddennum行inputnum列的二維矩陣
net.lw{2,1}=reshape(w2,outputnum,hiddennum);%更改矩陣的儲存格式
net.b{1}=reshape(B1,hiddennum,1);%1行hiddennum列，為隐含層的神經元門檻值
net.b{2}=reshape(B2,outputnum,1);
%網絡訓練
net=train(net,inputn,outputn);
an=sim(net,inputn);
error=sum(abs(an-outputn));%染色體對應的适應度值
% [m n]=size(inputn);
% A1=tansig(net.iw{1,1}.*inputn+repmat(net.b{1},1,n)); %需與main函數中激活函數相同
% A2=purelin(net.iw{2,1}.*A1+repmat(net.b{2},1,n)); %需與main函數中激活函數相同
% error=sumsqr(outputn-A2);
到這裡，我們已經得到了種群的初始化資訊，每個基因已經指派，并且每個個體的适應度值已經得到（适應度值是誤差，是以越小越好），那我們再通過種群的初始化資訊計算出最好的适應度值，找到那個個體。
%找最好的染色體
[bestfitness, bestindex]=min(individuals.fitness);%bestindex是最小值的索引（位置/某個個體），bestfitness的值為最小适應度值
bestchrom=individuals.chrom(bestindex,:); %最好的染色體，從10個個體中挑選到的
avgfitness=sum(individuals.fitness)/sizepop; %染色體的平均适應度(所有個體适應度和 / 個體數)
% 記錄每一代進化中最好的适應度和平均适應度
trace=[avgfitness bestfitness]; %trace矩陣，1行2列，avgfitness和bestfitness僅僅是數值
以上工作做完，我們就已經能夠把握這個種群的基本資訊了，但是肯定，她是第一代，不是最好的，那麼我們就要讓種群進行選擇，選出最好的個體，進行繁殖生子（交叉、變異），不斷地得到新的個體，再不斷地選擇……直到種群能得到最好的個體，才能最有效病準确地解決我們的問題，我們設定讓種群進化20代。
在這之前，我們要先來了解一下，如何讓代碼實作自然界中的“選擇”，并不是說個體的适應度值不好，就永遠不會被選擇到，真實情況應該是：适應度值好的個體，被選中的機率就會高一些，但是它還是有可能不會被選到的，這才叫公平。
我們已經把是個個體的适應度值随機初始化了，那麼我們可以把這個适應度值都轉成一個個機率，機率大的被選中機率就大一點，所有個體被選中的機率加起來就等于1。把十個個體集中放到一個輪盤裡，機率的大小對應着扇形面積的大小，然後我們産生0~1的随機數，随機數在哪個扇形區，哪個對應個體就會被選中。看代碼實作：
function ret=Select(individuals,sizepop)
% 本函數對每一代種群中的染色體進行選擇，以進行後面的交叉和變異
% individuals input : 種群資訊
% sizepop input : 種群規模
% ret output : 經過選擇後的種群
% 選擇操作即為決定整個遺傳算法種群走向的關鍵性步驟，遺傳算法的選擇操作
% 通常以較大的機率保留父代優良個體到子代中，以此使得種群往适應度值最佳的方向進化，
% 通俗的講就是：有若幹個備選方案，而且每個方案都有自己的潛力分值，但是在選擇的時候并不完全按照分值的高低來選
% 而是有一定的機率接受，分值高的接受機率高，分值較低接受的機率也低。
%根據個體适應度值進行排序
fitness1=10./individuals.fitness;%10為系數，取倒數得到新的數組fitness1，其值越大個體被選中的機率就越大
%.*或者./表示對矩陣裡每個元素都執行同樣操作
sumfitness=sum(fitness1);
sumf=fitness1./sumfitness;%每個個體遺傳到下一代的機率，1行sizepop列
index=[]; %僅僅是一個角标，當做記号即可
%采用輪盤賭法選擇新個體
for i=1:sizepop %轉sizepop次輪盤，%共選出sizepop個個體保留至下代
pick=rand; %産生随機數賦給pick
while pick==0 %若随機數是0，則重新産生随機數
pick=rand;
end
for j=1:sizepop %循環每個個體，找出此次的rank落在哪個個體的區間
pick=pick-sumf(j); %尋找rank落在哪個個體所在的區間，個體被選中的機率越大，其所占區間越大
if pick<0 %說明個體j的機率比較大（适應度值較小）
index=[index j]; %1行sizepop列的矩陣，存放被選中的個體标記
break; %尋找落入的區間，此次轉輪盤選中了染色體i，注意：在轉sizepop次輪盤的過程中，有可能會重複選擇某些染色體
end
end
end
%新的種群
individuals.chrom=individuals.chrom(index,:);
individuals.fitness=individuals.fitness(index);
ret=individuals;
如果算法已經能夠實作自然界中的“公平選擇”，那麼我們就讓它進化吧，不斷産生子代、變異，然後再一輪一輪地選擇。
%% 疊代求解最佳初始閥值和權值
% 進化開始
for i=1:maxgen
% 選擇
individuals=Select(individuals,sizepop);
% avgfitness=sum(individuals.fitness)/sizepop;%種群的平均适應度值
%交叉
individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop);
% 變異
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
以上代碼就已經實作了種群的進化和選擇，把不好的個體淘汰掉，把優秀的個體留下，下面來看交叉函數與變異函數：
function ret=Cross(pcross,lenchrom,chrom,sizepop)
%本函數完成交叉操作
% pcorss input : 交叉機率
% lenchrom input : 染色體的長度
% chrom input : 染色體群
% sizepop input : 種群規模
% ret output : 交叉後的染色體
for i=1:sizepop %每一輪for循環中，可能會進行一次交叉操作，染色體是随機選擇的，交叉位置也是随機選擇的，%但該輪for循環中是否進行交叉操作則由交叉機率決定（continue控制）
% 随機選擇兩個染色體進行交叉
pick=rand(1,2); %産生1行2列的随機數矩陣
while prod(pick)==0 %prod(A)檢測矩陣是否有0，有則重新生成随機數
pick=rand(1,2);
end
index=ceil(pick.*sizepop); %ceil(z)函數将輸入z中的元素取整,2.1則取3
% 交叉機率決定是否進行交叉
pick=rand;
while pick==0
pick=rand;
end
if pick>pcross
continue; %如果随機數比交叉機率大，則重新回到for語句，不執行下面的語句
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);%随機選到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(chrom(index(1),:)); %檢驗染色體1的可行性
flag2=test(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);%ceil(z)函數将輸入z中的元素取整,2.1則取3，得到随機一條染色體
% 變異機率決定該輪循環是否進行變異
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(index,pos)=chrom(index,pos)+(bound(pos,2)-chrom(index,pos))*fg;
else
chrom(index,pos)=chrom(index,pos)-(chrom(index,pos)-bound(pos,1))*fg;
end %變異結束
flag=test(chrom(index,:)); %檢驗染色體的可行性
end
end
ret=chrom;
好了，到這裡，我們已經通過代碼讓這個種群進化了20代，挑選出了對于有利于解決我們問題的個體，形象地對應起來就是，現在我們每個基因都有了最好的指派，當把它們放到網絡中之後，它們就能夠來解決問題了（預測作用）
%% 把最優初始閥值權值賦予網絡預測
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}=reshape(B2,outputnum,1);
回想一下此時此刻的網絡，已經有了權值門檻值，輸入層、隐含層和輸出層已經有明确的數值了，我們接下來再訓練好這個網絡就行了。
%% 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;
接着，再把改列印出來的圖形設定一下：
%% 遺傳算法結果分析
figure(1)
[r, c]=size(trace);
plot([1:r]',trace(:,2),'b--');
title(['适應度曲線 ' '終止代數＝' num2str(maxgen)]);
xlabel('進化代數');ylabel('适應度');
legend('平均适應度','最佳适應度');
disp('适應度 變量');
figure(2)
plot((output_test-test_simu)./test_simu,'-*');
title('神經網絡預測誤差百分比')
figure(3)
plot(error,'-*')
title('BP網絡預測誤差','fontsize',12)
ylabel('誤差','fontsize',12)
xlabel('樣本','fontsize',12)
完工，到這裡就能夠實作用遺傳算法來優化BP神經網絡了，最後運作得到圖形就能發現，加了遺傳算法之後網絡的預測能力就更強，誤差更小了。
如果要用代碼的朋友，不要複制上面的，因為我已經把它拆開了，應該有好幾個.m檔案的，最後，包括使用到的資料，我已經打包好了，代碼如果我有了解錯誤，歡迎朋友們糾正，謝謝！
工程下載下傳連結
（百度雲連接配接，密碼：gz73）
![深度學習與圍棋：為AlphaGo訓練深度神經網絡13.1.1 AlphaGo的網絡架構13.1.2 AlphaGo棋盤編碼器13.1.3 訓練AlphaGo風格的政策網絡[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)