通俗解釋matlab之遺傳演算法程式部分(二)
(1)程式怎麼開始
從哪裡開始程式比較好了？直接先主函式吧，然後再分著說：
%-------------函式說明----------------
% 主函式
%---------------------------------------
function main()
clear
clc
popsize = 100; %種群大小
chromlength = 10; %二進位制編碼長度
pc = 0.6; %交叉概率
pm = 0.001; %變異概率
pop = initpop(popsize,chromlength); %初始種群
for i=1:100
[objvalue] = cal_objvalue(pop); %計算適應度值（函式值）
fitvalue = objvalue;
[newpop] = selection(pop,fitvalue); %選擇操作
[newpop] = crossover(newpop,pc); %交叉操作
[newpop] = mutation(newpop,pm); %變異操作
pop = newpop; %更新種群
[bestindividual,bestfit]=best(pop,fitvalue);%尋找最優解
x2 = binary2decimal(bestindividual);
x1 = binary2decimal(newpop);
[y1] = cal_objvalue(newpop);
if mod(i,10)==0
figure;
fplot('10*sin(5*x)+7*abs(x-5)+10',[0 10]);
hold on;
title(['迭代次數為 n=' num2str(i)]);
plot(x1,y1,'*');
end
end
fprintf('the best X is --->>%5.2f\n'
fprintf('the best Y is --->>%5.2f\n',bestfit);
好了，這就是主函式，最後執行這個基本上什麼都有了，包括圖行什麼的都畫了，下面在分著說吧。。。
（2）關於二進位制種群怎麼生成
從上面的主程式可以看到，我們設定了100個個體（也就是100個初始化x值），二進位制編碼長度為10位，那麼先開始自然是生成這100個隨機個體了，注意是隨機的，不能使相同的，也就是說生成100個不同的二進位制編碼組合，像11000 10101這樣的，生成100組。怎麼辦，直接上程式：
%-------------函式說明----------------
% 初始化種群大小
% 輸入變數：
% popsize:種群大小
% chromlength：染色體長度--》轉化的二進位制長度
% 輸出變數：
% pop：種群
%---------------------------------------
function pop = initpop(popsize,chromlength)
pop = round(rand(popsize,chromlength));
很簡單，一句話搞定，說一下，關於rand（m,n）用法，就是生成m行n列的0~1之間隨機數，比如rand(3,4)為：
>> rand(3,4)
ans =
0.6028 0.1174 0.4242 0.2625
0.7112 0.2967 0.5079 0.8010
0.2217 0.3188 0.0855 0.0292
round什麼意思？四捨五入，簡單，這樣上面就變成了什麼了，round（ans）：
>> round(ans)
ans =
1 0 0 0
1 0 1 1
0 0 0 0
這樣當popsize=100，chromlength=10時就生成了對應的種群個體了。
（3）如何在把二進位制返回對應的十進位制
知道了二進位制怎麼返回對應範圍內的x值呢？前面講了演算法怎麼構造了，特別要注意的是二進位制的位數以及轉化完後的x值範圍，程式如下：
%-------------函式說明----------------
% 二進位制轉化十進位制函式
% 輸入變數：
% 二進位制種群
% 輸出變數：
% 十進位制數值
%---------------------------------------
function pop2 = binary2decimal(pop)
[px,py]=size(pop);
for i=1:py
pop1(:,i) = 2.^(py-i).*pop(:,i);
end
%sum(.,2)對行求和，得到列的向量
temp = sum(pop1,2);
pop2 = temp*10/1023;
輸入的為100組0、1編碼的二進位制，輸出的是x值。開始取一下種群大小size(pop)，顯然這裡py=10了，接著對每一位求和，就是pop1(:,i) = 2.^(py-i).*pop(:,i);這裡省略用了冒號，什麼意思了，就是對所有的行都是這個操作，冒號意思就是從1到100了。那麼就其中某個個體比如第一個吧，假設為11001 10010，那麼先進行這麼一個操作後就是什麼了？是不是就是對應為的0或1乘以2的對應次冪了，若果是1就管用，是0就不管用，那麼這個值就是：（2^0）*1+（2^1）*1+0+0+(2^4)*1+...,這樣就算出了一個值了，因為是10位編碼，所以這個數介於0~2^9即0~1023。那麼最大為多少？1023吧。temp = sum(pop1,2);對行求和吧，2表示對行，1表示對列，像
ans =
1 0 0 0
1 0 1 1
0 0 0 0
則>> sum(ans,1)
ans =
2 0 1 1
而>> sum(ans,2)
ans =
3
2
2
明白了吧，這樣temp就變成了有100行1列的值，且每個值都在0~1023之間變化，最後一行不解釋了，就是把它轉化為100組0~10之間的數值了。
（4）計算適應度函式值
這裡也就是根據100組x計算對應的100組y的值，簡單，根據上面的x值帶入就可以了：
%-------------函式說明----------------
% 計算函式目標值
% 輸入變數：
% 二進位制數值
% 輸出變數：
% 目標函式值
%---------------------------------------
function [objvalue]=cal_objvalue(pop)
x = binary2decimal(pop);
%轉化二進位制數為x變數的變化域範圍的數值
objvalue = 10*sin(5*x)+7*abs(x-5)+10;
（5）如何選擇新的個體
上面所有個體的函式值都計算出來了，存在objvalue中，此時它是不是也是100組y值呀，恩。那麼對於現有的隨機生成的100組x，怎麼來在選擇100組新的更好的x呢？這裡我們把選擇放在了交叉與變異之間了，都可以。如何選擇，就要構造概率的那個輪盤了，誰的概率大，是不是選擇的個體就會多一些？也就是現在的選擇就是100中選100個，最後出現的結果就是以前的100箇中最優的x有一個的話，選擇完後，可能就變成了5個這個x了，多餘的4個是不是相當於頂替了以前的不好的4個x值，這樣才能達到x總數100不變呀。
%-------------函式說明----------------
% 輸入變數：
% pop ： 二進位制種群
% fitvalue ： 適應度值
% 輸出變數：
% newpop：選擇以後的二進位制種群
%---------------------------------------
function [newpop] = selection(pop,fitvalue)
%構造輪盤
[px,py]=size(pop);
totalfit = sum(fitvalue);
p_fitvalue = fitvalue/totalfit;
p_fitvalue = cumsum(p_fitvalue);%概率求和排序
%-------
ms = sort(rand(px,1));%從小到大排列
fitin = 1;
newin = 1;
while newin<=px
if (ms(newin))<p_fitvalue(fitin)
newpop(newin,:)=pop(fitin,:);
newin=newin+1;
else fitin=fitin+1;
end
end
這一部分可能不太好理解。自己好好想想吧。前三句，求出每個個體被選擇的概率吧，第四句，求和排序是幹什麼的呢？比如現在假設一個概率可能是：
p_fitvalue =
0.0816
0.0408
0.1020
0.1224
0.0408
0.1429
0.1837
0.0204
0.0816
0.1837
那麼執行這一句p_fitvalue = cumsum(p_fitvalue)後就變為：
p_fitvalue =
0.0816
0.1224
0.2245
0.3469
0.3878
0.5306
0.7143
0.7347
0.8163
1.0000
也就是後面的值加到前面去在替換，這樣做有什麼意義呢？看後面吧。ms = sort(rand(px,1))這一句，生成100個0~1的隨機概率數，然後在從小到大排序，100次，好理解，輪盤賭不得賭100次嗎？為什麼排序呢？這和前面的cumsum求和了有關。好了真正的選擇開始了，while裡面，先對小的概率進行選擇，
if (ms(newin))<p_fitvalue(fitin)，就選擇一次，個體就多了一個。符合，newin概率的那個排序往上走，再符合，就再選擇一次，直到不符合了，fitin就往上走，這樣p_fitvalue(fitin)這部分是不是變大了一點，變大了再看看這個值能被選擇幾次。其實每一次迴圈就是對於某一個個體來說，看能選擇幾次。比如假設fitin=9時，此時假設p_fitvalue(9)=0.5，被選擇完了，下一次時fitin=10，我們假設離譜點，p_fitvalue(10)=0.7的話，那麼在下一次去ms（newin）的話(假設現在已經選擇了newin=30個個體了)，是不是一定為選擇呀，然後第31個個體就是p_fitvalue(10)這個個體了，在32個，因為ms（newin）也是從小到大排序的，所以這次選擇只是增加了一點點吧，假設為0.52（肯定大於0.5吧，不然也不會進來），0.52<0.7，又選擇p_fitvalue(10)這個個體吧，這樣一直下去，33,34，，，好，來到了第50個個體假設。此時ms（50）=0.71了，不行了，是不是不行了，fitin是不是加1變為11了，如果p_fitvalue(11)較p_fitvalue(10)變化不大，比如只變為了0.72，那麼第11個個體的被選擇的概率是不是隻是0.02呀，這麼小點，也只能是就選擇那麼2回就下一個了吧，而我們看看第10個個體，從0.5變到了0.7，它的概率是不是0.2呀，哇，0.2呀，好大呀，難怪能從第32個個體一直選擇到第50個個體呀。好了到這裡應該能明白了吧，再不明白就沒招了，自己好好想吧。那麼這麼一系列迴圈後，是不是100個個體就被選擇出來了呀，想想選擇的結果是什麼樣子的呢？100個個體是不是會出現為：第2個第2個，第3個第3個第3個，第5個第5個第5個第5個第5個，第6個，第10個第10個第10個，等等等，總共100個個體吧（那些重複的數是在以前種群中的編號）。好多呀，應該講明白了選擇這塊。
（6）怎麼交叉
選擇完後是不是就有新的相對優秀的100個個體了，那麼對它們我們在來進行交叉變異才能形成新的個體，畢竟這100個個體還是先前100個個體裡面挑好的挑出來的吧。
%-------------函式說明----------------
% 輸入變數：
% pop:二進位制的父代種群數
% pc：交叉的概率
% 輸出變數：
% newpop：交叉後的種群數
%---------------------------------------
function [newpop]=crossover(pop,pc)
[px,py]=size(pop);
newpop = ones(size(pop));
for i=1:2:px-1
if(rand<pc)
cpoint = round(rand*py);
newpop(i,:) = [pop(i,1:cpoint),pop(i+1,cpoint+1:py)];
newpop(i+1,:) = [pop(i+1,1:cpoint),pop(i,cpoint+1:py)];
else
newpop(i,:)=pop(i,:);
newpop(i+1,:)=pop(i+1,:);
end
end
這裡涉及到了一個概率，就是交叉的概率pc，因為交叉不是一定發生，所以用一個概率來限制，部分還是要保持原來不變的。程式開始先生產一個與原種群同大小的全1矩陣newpop，然後迴圈for i=1:2:px-1，想想為什麼間距是2呢？簡單，每兩個交叉吧，就像自然界偉大的規律一樣，當然也有例外，不考慮吧。然後生產隨機0~1數rand，在一比較就知道是不是該交叉，這裡的pc = 0.6，那麼想想每次rand的話，由於隨機性，是不是它產生小於0.6的概率就是0.6呀。好了選擇了交叉，就執行下面操作，其中還有個交叉點cpoint 計算，能明白吧，不解釋，後面兩句也能看懂。如果概率比0.6大的話（也就是0.4的概率rand出了0.6~1之間的數），就不選擇交叉，直接賦值就行了吧。這樣知道結束，有個問題，萬一px為奇數怎麼辦？是不是相當於有一個光棍呀，是呀，咋辦的呢？我也不知道，但是我們這裡為100，還好是剛好1配1呀，很人性，大家程式設計序也要人性化一點呀，雖然那是數字，沒有生命，但是體現了你的那個什麼，道德修養，對世界的愛吧~_~。
（7）關於變異
好了交叉完後就生成了新一代的種群了吧，那麼在對這個種群再一次進行變異，實際情況下交叉變異應該是同時發生的吧，這裡先後無所謂了。那麼怎麼變異呢？
%-------------函式說明----------------
% 輸入變數：
% pop ： 二進位制種群
% pm ： 變異概率
% 輸出變數：
% newpop ： 變異以後的種群
%---------------------------------------
function [newpop] = mutation(pop,pm)
[px,py] = size(pop);
newpop = ones(size(pop));
for i=1:px
if(rand<pm)
mpoint = round(rand*py);
if mpoint<=0
mpoint=1;
end
newpop(i,:) = pop(i,:);
if newpop(i,mpoint)==0
newpop(i,mpoint)=1;
else newpop(i,mpoint)=0;
end
else
newpop(i,:)=pop(i,:);
end
end
個人感覺變異更簡單，就是找個點（高階一點的找多個點），把它的值（0或1）取反不就得了嘛，程式就像上面的那樣，講到這裡了感覺真不需要解釋了，略過。
（8）最後一點了
至此，主體部分算是完了吧，剩下的就是顯示、處理部分了，回到主程式去看吧，接下來就是更新種群吧，然後再變換十進位制去，提取最好的y值和對應的x值，這部分還有個求最優適應度的函式：
%-------------函式說明----------------
% 輸入變數：
% pop : 種群
% fitvalue ： 種群適應度
%
% 輸出變數：
% bestindividual ： 最佳個體（二進位制個體）
% bestfit ： 最佳適應度值
%---------------------------------------
function [bestindividual,bestfit]=best(pop,fitvalue)
[px,py]=size(pop);
bestindividual = pop(1,:);
bestfit = fitvalue(1);
for i=2:px
if fitvalue(i)>bestfit
bestindividual = pop(i,:);
bestfit = fitvalue(i);
end
end
感覺也沒有什麼號解釋的，就是比較唄，最優值在bestfit裡面，如果有比它大的，就更新，否則略過，不解釋。
最後是在matlab環境下畫出來，if mod(i,10)==0求餘操作吧，也就是每迭代10代我畫一次圖，觀察結果吧，後面的fprintf，打印出來，很簡單，自己看吧。。。
至此，整個遺傳演算法的大致原理講解完了，程式的各部分程式碼都在這裡了，把每個小部分都編寫一個m檔案存起來再一起執行就ok了。所有的程式排列以及這部分實驗結果在下一個裡面再貼出來吧~_~