一、簡介

1 遺傳演算法概述
遺傳演算法（Genetic Algorithm，GA）是進化計算的一部分，是模擬達爾文的遺傳選擇和自然淘汰的生物進化過程的計算模型，是一種通過模擬自然進化過程搜尋最優解的方法。該演算法簡單、通用，魯棒性強，適於並行處理。

2 遺傳演算法的特點和應用
遺傳演算法是一類可用於複雜系統優化的具有魯棒性的搜尋演算法，與傳統的優化演算法相比，具有以下特點：
（1）以決策變數的編碼作為運算物件。傳統的優化演算法往往直接利用決策變數的實際值本身來進行優化計算，但遺傳演算法是使用決策變數的某種形式的編碼作為運算物件。這種對決策變數的編碼處理方式，使得我們在優化計算中可借鑑生物學中染色體和基因等概念，可以模仿自然界中生物的遺傳和進化激勵，也可以很方便地應用遺傳操作運算元。
（2）直接以適應度作為搜尋資訊。傳統的優化演算法不僅需要利用目標函式值，而且搜尋過程往往受目標函式的連續性約束，有可能還需要滿足“目標函式的導數必須存在”的要求以確定搜尋方向。遺傳演算法僅使用由目標函式值變換來的適應度函式值就可確定進一步的搜尋範圍，無需目標函式的導數值等其他輔助資訊。直接利用目標函式值或個體適應度值也可以將搜尋範圍集中到適應度較高部分的搜尋空間中，從而提高搜尋效率。
（3）使用多個點的搜尋資訊，具有隱含並行性。傳統的優化演算法往往是從解空間的一個初始點開始最優解的迭代搜尋過程。單個點所提供的搜尋資訊不多，所以搜尋效率不高，還有可能陷入區域性最優解而停滯；遺傳演算法從由很多個體組成的初始種群開始最優解的搜尋過程，而不是從單個個體開始搜尋。對初始群體進行的、選擇、交叉、變異等運算，產生出新一代群體，其中包括了許多群體資訊。這些資訊可以避免搜尋一些不必要的點，從而避免陷入區域性最優，逐步逼近全域性最優解。
（4） 使用概率搜尋而非確定性規則。傳統的優化演算法往往使用確定性的搜尋方法，一個搜尋點到另一個搜尋點的轉移有確定的轉移方向和轉移關係，這種確定性可能使得搜尋達不到最優店，限制了演算法的應用範圍。遺傳演算法是一種自適應搜尋技術，其選擇、交叉、變異等運算都是以一種概率方式進行的，增加了搜尋過程的靈活性，而且能以較大概率收斂於最優解，具有較好的全域性優化求解能力。但，交叉概率、變異概率等引數也會影響演算法的搜尋結果和搜尋效率，所以如何選擇遺傳演算法的引數在其應用中是一個比較重要的問題。
綜上，由於遺傳演算法的整體搜尋策略和優化搜尋方式在計算時不依賴於梯度資訊或其他輔助知識，只需要求解影響搜尋方向的目標函式和相應的適應度函式，所以遺傳演算法提供了一種求解複雜系統問題的通用框架。它不依賴於問題的具體領域，對問題的種類有很強的魯棒性，所以廣泛應用於各種領域，包括：函式優化、組合優化生產排程問題、自動控制
、機器人學、影象處理（影象恢復、影象邊緣特徵提取......)、人工生命、遺傳程式設計、機器學習。

3 遺傳演算法的基本流程及實現技術
基本遺傳演算法（Simple Genetic Algorithms,SGA）只使用選擇運算元、交叉運算元和變異運算元這三種遺傳運算元，進化過程簡單，是其他遺傳演算法的基礎。

3.1 遺傳演算法的基本流程
通過隨機方式產生若干由確定長度（長度與待求解問題的精度有關）編碼的初始群體；
通過適應度函式對每個個體進行評價，選擇適應度值高的個體參與遺傳操作，適應度低的個體被淘汰；
經遺傳操作（複製、交叉、變異）的個體集合形成新一代種群，直到滿足停止準則（進化代數GEN>=?）；
將後代中變現最好的個體作為遺傳演算法的執行結果。

其中，GEN是當前代數；M是種群規模，i代表種群數量。

3.2 遺傳演算法的實現技術
基本遺傳演算法（SGA）由編碼、適應度函式、遺傳運算元（選擇、交叉、變異）及執行引數組成。
3.2.1 編碼
（1）二進位制編碼
二進位制編碼的字串長度與問題所求解的精度有關。需要保證所求解空間內的每一個個體都可以被編碼。
優點：編、解碼操作簡單，遺傳、交叉便於實現
缺點：長度大
（2）其他編碼方法
格雷碼、浮點數編碼、符號編碼、多引數編碼等
3.2.2 適應度函式
適應度函式要有效反映每一個染色體與問題的最優解染色體之間的差距。
3.2.3選擇運算元

3.2.4 交叉運算元
交叉運算是指對兩個相互配對的染色體按某種方式相互交換其部分基因，從而形成兩個新的個體；交叉運算是遺傳演算法區別於其他進化演算法的重要特徵，是產生新個體的主要方法。在交叉之前需要將群體中的個體進行配對，一般採取隨機配對原則。
常用的交叉方式：
單點交叉
雙點交叉（多點交叉，交叉點數越多，個體的結構被破壞的可能性越大，一般不採用多點交叉的方式）
均勻交叉
算術交叉
3.2.5 變異運算元
遺傳演算法中的變異運算是指將個體染色體編碼串中的某些基因座上的基因值用該基因座的其他等位基因來替換，從而形成一個新的個體。

就遺傳演算法運算過程中產生新個體的能力方面來說，交叉運算是產生新個體的主要方法，它決定了遺傳演算法的全域性搜尋能力；而變異運算只是產生新個體的輔助方法，但也是必不可少的一個運算步驟，它決定了遺傳演算法的區域性搜尋能力。交叉運算元與變異運算元的共同配合完成了其對搜尋空間的全域性搜尋和區域性搜尋，從而使遺傳演算法能以良好的搜尋效能完成最優化問題的尋優過程。

3.2.6 執行引數

4 遺傳演算法的基本原理
4.1 模式定理

4.2 積木塊假設
具有低階、定義長度短，且適應度值高於群體平均適應度值的模式稱為基因塊或積木塊。
積木塊假設：個體的基因塊通過選擇、交叉、變異等遺傳運算元的作用，能夠相互拼接在一起，形成適應度更高的個體編碼串。
積木塊假設說明了用遺傳演算法求解各類問題的基本思想，即通過積木塊直接相互拼接在一起能夠產生更好的解。

二、原始碼

%Run NAA;
function OUT = MO_NAA( funcName_adjustInd, funcName_fitness, userObj, controlParams)     
%user parameters;
NOBJ = userObj.NOBJ;
popSize = userObj.popSize;
dimension = userObj.dimension;
bounds = userObj.bounds;
types = userObj.types;
MAXGEN = userObj.MAXGEN;
MAXFUNEVALS = userObj.MAXFUNEVALS;

CounterGEN = userObj.CounterGEN;
CounterFES = userObj.CounterFES;

%control parameters;
shelterNum = controlParams.shelterNum;
shelterCap = zeros(1, shelterNum) + controlParams.shelterCap;
scale_local = controlParams.scale_local;
Cr_local = controlParams.Cr_local;
alpha = controlParams.alpha;
Cr_global = controlParams.Cr_global;

bounceBack = controlParams.bounceBack;

%initiallize the population;
[swarm] = MO_initializeSwarm(popSize, dimension, bounds);  %initialize the individuals. Each individual represents a solution;  %% is a popSize*dimension matrix
%initialize the shelters;
[swarm, swarmFitnesses, shelterIndexes, shelterLeaders, normalizedLeaderFitnesses, userObj]  = MO_formShelters(swarm, dimension, bounds, types, popSize, shelterNum, shelterCap, funcName_adjustInd, funcName_fitness, userObj);
CounterFES = userObj.CounterFES;

PFront = zeros(popSize,NOBJ);
PSet = zeros(popSize,dimension);

%start the iteration process;
for i=1:MAXGEN   

   indNumbers_shelters = MO_calculateIndNumberInShelters(shelterNum, shelterIndexes);  %Calculate the number of individuals present in the shelters;
   
   %scan the individuals to do the search;
   for p=1:popSize
        shelterIndex = shelterIndexes(p);
        
        if(shelterIndex >0) %if the individual is currently in a shelter;
            theta = normalizedLeaderFitnesses(1,shelterIndex); %get the shelter quality value;  %% size(theta) = 1*1
            [swarm, shelterIndexes, swarmFitnesses, userObj, P, S] = MO_decideLeave(bounds, types, popSize, dimension, shelterIndex, theta, indNumbers_shelters(shelterIndex), shelterIndexes, swarm, p, shelterLeaders, swarmFitnesses, funcName_adjustInd, funcName_fitness, controlParams, userObj);
            CounterFES = userObj.CounterFES;
            
            PFront(p,:) = P;   %% pareto front
            PSet(p,:) = S;    %%  state variable for pareto front 
        
        else  %if the cockroach is currently doing the exploring;
            theta = zeros(1,1);
            r_c = 0;
            numberInCorner = 0;
            if (shelterNum>0)
                r_c = unidrnd(shelterNum);
                theta = normalizedLeaderFitnesses(1,r_c);  %% size(theta) = NOBJ*1
                numberInCorner = indNumbers_shelters(r_c);
            end
            
            [swarm, shelterIndexes, swarmFitnesses, userObj, P, S] = MO_decideEnter(bounds, types, popSize, dimension, theta, r_c, numberInCorner, shelterIndexes, swarm, p, shelterLeaders, swarmFitnesses, funcName_adjustInd, funcName_fitness, controlParams, userObj);
            CounterFES = userObj.CounterFES;

            PFront(p,:) = P;   %% pareto front
            PSet(p,:) = S;    %%  state variable for pareto front 
        end       
                      
   end
   function [newInd, newUserObj] = MO_adjustInd_test(ind, userObj)
NOBJ = userObj.NOBJ;
popSize = userObj.popSize;  
dimension = userObj.dimension;


k=6;
n=2*combntns(k,2);
t_number=3;
car_tnumber=100;
E_Demand=-2.0;
E_Time=-2.5;
C_mv=0.007;
C_mp=0.0013;
C_v=0.012;
p_charging=20;
P0=0.3;
B_v=30;
B_v_LT=0.5;
B_v_UT=0.95;
R_dc=0.1;
PEC=0.0032;
Z_k = [24   22  21	22   22  21];
a_kt = [5,6,3,2,4,7;7,6,2,1,6,4;1,4,5,5,6,4];
PEL=0.15; 
R_vdc=0.1;



X=ind;

P_kj=X(1,1:n); 
P_kjt(1,:)=X(1,n+1:2*n); 
P_kjt(2,:)=X(1,2*n+1:3*n); 
P_kjt(3,:)=X(1,3*n+1:4*n); 

%% calculated parameters

%% t0
t0_lower=10;
t0_upper=30;
t0_kj=t0_lower + (t0_upper-t0_lower).*rand([t_number n]);

%% t0ktv
t0ktv_lower=4;
t0ktv_upper=10;
t0_ktv=t0ktv_lower + (t0ktv_upper-t0ktv_lower).*rand([t_number k car_tnumber]);

%% D0
%different area
D01_lower=2;
D01_upper=5;
D02_lower=6;
D02_upper=9;
D03_lower=10;
D03_upper=16;
D0_kj1=round(D01_lower + (D01_upper-D01_lower).*rand([3 n/3]));
D0_kj2=round(D02_lower + (D02_upper-D02_lower).*rand([3 n/3]));
D0_kj3=round(D03_lower + (D03_upper-D03_lower).*rand([3 n/3]));
D0_kj=[D0_kj1,D0_kj2,D0_kj3];
%diferent time
d1=D0_kj(1,:)+2;
d2=D0_kj(2,:)+10;
d3=D0_kj(3,:)+4;
D0_kj=[d1;d2;d3];

%% soc0
soctv_lower=0.3;
soctv_upper=0.9;
SOC_tv=soctv_lower + (soctv_upper-soctv_lower).*rand([t_number car_tnumber]);


%% t_kj
for t=1:3
    for i=1:n
                 t_kj(t,i)=((((P_kjt(t,i)-P0)/P0)*E_Time*t0_kj(t,i))+t0_kj(t,i));
    end
end
for t=1:3
    for i=1:k
                 t_kjnew(t,i)=t_kj(t,i);
                 t0_kjnew(t,i)=t0_kj(t,i);
    end
end


%% D_kj
for t=1:3
    for i=1:n
                    D_kj(t,i)=((((P_kj(1,i)-P0)/P0)*E_Demand)*D0_kj(t,i)+D0_kj(t,i));
    end
end


%% t_ktv
for t=1:3
        for station=1:k
            for car=1:car_tnumber
                t_ktv(t,station,car)=t_kjnew(t,station)./t0_kjnew(t,station).*t0_ktv(t,station,car);   
            end
        end
end






%% check constraints.



%% constraint (23)-(25)
%(23)
for t = 1:3
    for j = 1:n
        if D_kj(t,j)>=D0_kj(t,j)+((D0_kj(t,j)*(P_kj(1,j)-P0)*E_Demand)/(P0))-0.5;
            success_23(t,j)=1;
        else
            success_23(t,j)=0;
        end
    end
end


% 
% 
%(24)
for t = 1:3
    for j = 1:n
        if D_kj(t,j)<=D0_kj(t,j)+((D0_kj(t,j)*(P_kj(1,j)-P0)*E_Demand)/(P0))+0.5;
            success_24(t,j)=1;
        else
            success_24(t,j)=0;
        end
    end
end
% 
% 
%(25)
for t = 1:3
    for j = 1:n
        if D0_kj(t,j)+((D0_kj(t,j)*(P_kj(1,j)-P0)*E_Demand)/(P0))>=0;
            success_25(t,j)=1;
        else
            success_25(t,j)=0;
        end
    end
end

三、執行結果

四、備註

版本：2014a