智能優化算法


Matlab

下標從1開始!!

  • ~=:相當於不等於!=

  • eps表示的是一個數可以分辨的最小精度,返回1.0和下一個精度最高的雙精度浮點數的差值, 即2^(-52)。

  • Inf和-Inf分別代表正無窮和負無窮

  • y = length(x) 函數計算指定向量或矩陣的長度y。如果參數變量x是向量,則返回其長度;如果參數變量是非空矩陣,則length(x)與max(size(x))等價

矩陣取值

  • v(m , : ):取第m行
  • v( : ,m):取第m列
  • v(x,y):取第x行第y列

矩陣計算

  • [SortFit1,Index]=sort(Fit1):對Fit進行排序,排序結果放入SortFit1矩陣,結果每位上的元素在原來列上的序號放入Index矩陣
  • [aa,bb]=size(A):aa=行數,bb=列數
  • cumsum(x):對矩陣x進行逐列累加,例:[a,b,c,d]=>[a,a+b,a+b+c,a+b+c+d]
  • sum(x):對矩陣x進行逐列求和,即每列之和

取整函數

  • fix朝零方向取整,如fix(-1.3)=-1; fix(1.3)=1;
  • floor,顧名思義,就是地板,所以是取比它小的整數,即朝下取整,如floor(-1.3)=-2; floor(1.3)=1;floor(-1.8)=-2,floor(1.8)=1
  • ceil,與floor相反,它的意思是天花板,也就是取比它大的最小整數,即朝上取整,如ceil(-1.3)=-1; ceil(1.3)=2;ceil(-1.8)=-1,ceil(1.8)=2
  • round四舍五入到最近的整數,如round(-1.3)=-1;round(-1.52)=-2;round(1.3)=1;round(1.52)=2

生成隨機數函數

  • rand 生成均勻分布的偽隨機數。分布在(0~1)之間
    • rand(m,n)生成m行n列的均勻分布的偽隨機數
  • randn 生成標准正態分布的偽隨機數(均值為0,方差為1)
  • randi 生成均勻分布的偽隨機整數
    • randi(iMax)在閉區間[1,iMax]生成均勻分布的偽隨機整數
    • randi(iMax,m,n)在開區間[1,iMax]生成mXn型隨機矩陣
    • r = randi([iMin,iMax],m,n)在開區間[iMin,iMax]生成m*n型隨機矩陣
  • randperm 生成整數的隨機排列

排序函數

B = sort(A)升序A 的元素進行排序。

  • 如果 A 是向量,則 sort(A) 對向量元素進行排序。
  • 如果 A 是矩陣,則 sort(A) 會將 A 的列視為向量並對每列進行排序,從上到下,從小到大
  • 如果 A 是多維數組,則 sort(A) 會沿大小不等於 1 的第一個數組維度計算,並將這些元素視為向量

復制函數

B = repmat(A,x,y)生成重復數矩陣

  • 若A是一個數,則生成x*y的矩陣,全是A
  • 若A是一個矩陣,將矩陣A復制2行3列

智能優化算法

代表人物匯總

算法名 英文名 代表人物
遺傳算法GA Genetic Algorithm J.H.Holand
差分進化算法DE Differential evolution Storn
免疫算法IA Immune algorithm Burnet
蟻群算法ACO Ant Colony optimization M.Dorigo,V.Maniezzo,A.Colorni
粒子群算法PSO Particle swarm optimization James Kennedy,Rusell Eberhart
模擬退火算法SA Simulated annealing Metropolis
禁忌搜索算法TS Tabu Search Glover
神經網絡算法NN Neural Newwork McCulloch,Pitts,J.J.Hopfield

算法特點匯總

算法名字 特點 優缺點
遺傳算法 群體搜索策略和簡單的遺傳算子 全局搜索能力強,局部搜索能力較弱,早熟,算法收斂性無法保證
差分進化算法
免疫算法
蟻群算法
粒子群算法 收斂速度快但容易陷入局部最優解
模擬退火算法
禁忌搜索算法
神經網絡算法

遺傳算法GA

​ 遺傳算法模擬生物在自然環境中的遺傳和進化的過程,從而形成自適應全局優化搜索算法,它借用了生物遺傳學的觀點,通過自然選擇,遺傳,變異等機制,實現各個個體適應性的提高

名詞解釋

遺傳學術語 遺傳算法術語
群體 可行解集
個體 可行解
染色體 可行解的編碼
基因 可行解編碼的分量
基因形式 遺傳編碼
適應度 評價函數值
選擇 選擇操作
交叉 交叉操作
變異 變異操作

關鍵參數

  • 群體規模Np:群體規模將影響遺傳優化的最終結果以及遺傳算法的執行效率。一般取10~200。
  • 交叉概率Pc:交叉概率控制着交叉操作被使用的頻度。一般取0.25~1.00。
  • 變異概率Pm:變異的主要目的是保持群體的多樣性,一般地頻度的變異可防止群體中重要基因的丟失。一般取0.001~0.1。
  • 遺傳運算的終止進化代數G:終止遺傳代數表示遺傳算法運行到指定的進化代數之后就停止運行。一般取100~1000。

主要流程

選擇

輪盤賭選擇(Roulette Wheel Selection)

一種回放式隨機采樣方法。每個個體進入下一代的概率等於它的適應度值與整個種群中個體適應度值和的比例

輪盤賭選擇法是最簡單也是最常用的選擇方法,在該方法中,各個個體的選擇概率和其適應度值成比例,適應度越大,選中概率也越大。但實際在進行輪盤賭選擇時個體的選擇往往不是依據個體的選擇概率,而是根據“累積概率”來進行選擇。

例子:

  1. 假設群體 A:a, b, c, d, e
  2. 適應度 F1:0.45, 0.1, 0.1, 0.15, 0.2
  3. 對F1進行累加操作得到F2(累積概率):0.45, 0.55, 0.65, 0.8, 1.0
  4. 生成排序隨機隊列MS:0.2, 0.3, 0.4 , 0.7 ,0.8 ,0.9
  5. 選擇完成群體 B:a, a, a, d, e, e

合並新舊種群取最優

該選擇策略,是在交配產生的新子代和父代中,通過合並子代和父代,並且按照個體的適應度進行排序,選擇適應度最佳的前 NP 個個體作為新的種群,達到更好的收斂的效果。

例子:

  1. 父代個體parent : a,b,c
  2. 父代適應度parent_obj : 1 , 5 , 3
  3. 子代個體son : d , e , f
  4. 子代適應度son_obj : 2 6 4
  5. 合並個體total = a,b,c,d,e,f
  6. 合並適應度total_obj : 1,5,3,2,6,4
  7. 適應度排序order : e,b,f,c,d,a
  8. 前 NP=3個作為下一代:save = e,b,f

交叉

單點交叉

單點交叉通過選取兩條染色體,在隨機選擇的位置點上進行分割並交換右側的部分,從而得到兩個不同的子染色體。

例子:

(空格處代表該處為隨機生成的交叉點,往后的基因需要進行交換操作)

待交叉染色體: 11111 22222

​ 33333 44444

交換后染色體: 11111 44444

​ 33333 22222

君主交叉

君主交叉就是在種群內選擇出最優個體,用最優個體作為君主染色體,並隨機生成交叉位點,然后將君主染色體的特定位點上的基因遺傳給子代染色體

例子:

(中間空格分隔的基因代表隨機生成的交叉位點,該位點上的基因需要遺傳給子代染色體)

待交叉染色體: 101000 1 101

​ 011001 0 011

交換后染色體: 101000 1 101

​ 011001 1 011

多點奇偶交叉

  • 先規定一個交叉率,然后循環奇數位
  • 每次生成一個隨機數,判斷是否小於交叉率
    • 如果是,則進行交叉
      • 對二進制位上每一位取隨機數,如果等於1,則與偶數位同位置的數進行交換
    • 如果不是,則跳過交叉

變異

通過變異率計算第i個個體的j個基因二進制位上會發生變異,即取反

域內隨機值

描述:

對新種群內的所有染色體每個基因進行變異操作,先通過隨機值和變異率去判斷該基因是否會發生變異,再通過隨機值在基因的取值范圍內進行改變。

例子:

  1. 在所有染色體中隨機挑選 NP * Pm 個染色體進行變異
  2. 假設挑選出染色體h1:11000011
  3. 隨機挑選 L(染色體長度,基因個數) * Pm 個基因進行變異
  4. 挑選下標:1,2,4
  5. h2:00101011(下標1,2,4處基因取反)

代碼

奇偶交叉GA1

%%%%%%%%%%%%%%%%%%%%標准遺傳算法求函數極值%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化參數%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %清除所有變量
close all;              %清圖
clc;                    %清屏
NP=50;                  %種群數量
L=20;                   %二進制數串長度
Pc=0.8;                 %交叉率
Pm=0.1;                 %變異率
G=100;                  %最大遺傳代數
Xs=10;                  %上限
Xx=0;                   %下限
f=randi([0,1],NP,L);    
%%%%%%%%%%%%%%%%%%%%%%%%%遺傳算法循環%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:G            
    %%%%%%%%%%%%將二進制解碼為定義域范圍內十進制%%%%%%%%%%%%%%
    for i=1:NP         
        U=f(i,:);       
        m=0;
        for j=1:L      
            m=U(j)*2^(j-1)+m;
        end
        x(i)=Xx+m*(Xs-Xx)/(2^L-1); 
        Fit(i)= func1(x(i));
    end       
    maxFit = max(Fit);           %最大值
    minFit = min(Fit);           %最小值
    rr = find(Fit==maxFit);
    fBest = f(rr(1,1),:);        %歷代最優個體   
    xBest = x(rr(1,1));
    Fit = (Fit-minFit)/(maxFit-minFit);	%歸一化適應度值
    %%%%%%%%%%%%%%%%%%基於輪盤賭的復制操作%%%%%%%%%%%%%%%%%%%
    sum_Fit=sum(Fit);					%逐列求和
    fitvalue=Fit./sum_Fit;				%每個元素計算適應度百分比
    fitvalue=cumsum(fitvalue);			%逐列累加適應度百分比,求得累積概率
    ms=sort(rand(NP,1));				%生成排序隨機隊列
    fiti=1;
    newi=1;
    while newi<=NP						%選擇完成新群體,
        if (ms(newi))<fitvalue(fiti)	%當隨機隊列概率小於當前累積概率,則復制該個體為新群體上的新個體,新群體個體序號+1
            nf(newi,:)=f(fiti,:);
            newi=newi+1;
        else							%當隨機隊列概率大於等於當前累積概率,則累積概率數組序號+1
            fiti=fiti+1;
        end
    end   
    %%%%%%%%%%%%%%%%%%%%%%基於概率的交叉操作%%%%%%%%%%%%%%%%%%
    for i=1:2:NP						%1,3,5...奇數位遍歷
        p=rand;
        if p<Pc							%如果小於交叉率
            q=randi([0,1],1,L);			%生成1到L的數組,每位上隨機取0或1
            for j=1:L			
                if q(j)==1;				%若該位上是1,則與i+1上相同位上數進行swap
                    temp=nf(i+1,j);
                    nf(i+1,j)=nf(i,j);
                    nf(i,j)=temp;
                end
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%基於概率的變異操作%%%%%%%%%%%%%%%%%%%%%%%
    i=1;
    while i<=round(NP*Pm)		  	%迭代數=變異率x總數=變異染色體數量
        h=randi([1,NP],1,1);      	%隨機選取一個需要變異的染色體
        for j=1:round(L*Pm)         %在一個染色體上隨機變異一定數量的基因
            g=randi([1,L],1,1);   	%計算隨機需要變異的基因在第幾位
            nf(h,g)=~nf(h,g);		%將其取反=即是發生變異
        end
        i=i+1;
    end
    f=nf;
    f(1,:)=fBest;                   %保留最優個體在新種群中
    trace(k)=maxFit;                %歷代最優適應度
end
xBest;                              %最優個體

君主交叉GA2

%%%%%%%%%%%%%%%%%%%%實值遺傳算法求函數極值%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                           %清除所有變量
close all;                           %清圖
clc;                                 %清屏
D=10;                                %基因數目    
NP=100;                              %染色體數目
Xs=20;                               %上限          
Xx=-20;                              %下限
G=100;                               %最大遺傳代數
f=zeros(D,NP);                       %初始種群賦空間
nf=zeros(D,NP);                      %子種群賦空間
Pc=0.8;                              %交叉概率
Pm=0.1;                              %變異概率

f=rand(D,NP)*(Xs-Xx)+Xx;             %隨機獲得初始種群
%%%%%%%%%%%%%%%%%%%%%%按適應度升序排列%%%%%%%%%%%%%%%%%%%%%%%
for np=1:NP
    FIT(np)=func2(f(:,np));
end
[SortFIT,Index]=sort(FIT);                            
Sortf=f(:,Index);
%%%%%%%%%%%%%%%%%%%%%%%遺傳算法循環%%%%%%%%%%%%%%%%%%%%%%%%%%
for gen=1:G
    %%%%%%%%%%%%%%采用君主方案進行選擇交叉操作%%%%%%%%%%%%%%%%
    Emper=Sortf(:,1);                      	%君主染色體
    NoPoint=round(D*Pc);                   	%每次交叉基因的個數=round(基因數x交叉率)
    PoPoint=randi([1 D],NoPoint,NP/2);     	%交叉基因的位置,NoPoint x NP/2的隨機矩陣
    nf=Sortf;
    for i=1:NP/2
        nf(:,2*i-1)=Emper;					%奇數位=君主
        nf(:,2*i)=Sortf(:,2*i);				%偶數位=原基因	
        for k=1:NoPoint
            nf(PoPoint(k,i),2*i-1)=nf(PoPoint(k,i),2*i);
            nf(PoPoint(k,i),2*i)=Emper(PoPoint(k,i));
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%變異操作%%%%%%%%%%%%%%%%%%%%%%%%%
    for m=1:NP
        for n=1:D
            r=rand(1,1);	%隨機生成1個數在(0,1]之間
            if r<Pm			%若隨機概率小於變異率則發生變異
                nf(n,m)=rand(1,1)*(Xs-Xx)+Xx;	%即用概率乘以區間差+下限重新生成新的值
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%子種群按適應度升序排列%%%%%%%%%%%%%%%%%%
    for np=1:NP 
          NFIT(np)=func2(nf(:,np));   	%計算子種群適應度
    end
    [NSortFIT,Index]=sort(NFIT);        %子種群排序   
    NSortf=nf(:,Index);					%子種群
    %%%%%%%%%%%%%%%%%%%%%%%%%產生新種群%%%%%%%%%%%%%%%%%%%%%%%%%%
    f1=[Sortf,NSortf];                	%子代和父代合並
    FIT1=[SortFIT,NSortFIT];       		%子代和父代的適應度值合並
    [SortFIT1,Index]=sort(FIT1);    	%適應度按升序排列
    Sortf1=f1(:,Index);               	%按適應度排列個體
    SortFIT=SortFIT1(1:NP);         	%取前NP個適應度值
    Sortf=Sortf1(:,1:NP);             	%取前NP個個體
    trace(gen)=SortFIT(1);           	%歷代最優適應度值
end
Bestf=Sortf(:,1);                     	%最優個體 
trace(end)                            	%最優值

GA算法求解旅行商問題

%%%%%%%%%%%%%%%%%%%%%%%%%遺傳算法解決TSP問題%%%%%%%%%%%%%%%%%%%%%%%
clear all;                      %清除所有變量
close all;                      %清圖
clc;                            %清屏
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
    3238 1229;4196 1044;4312  790;4386  570;3007 1970;2562 1756;...
    2788 1491;2381 1676;1332  695;3715 1678;3918 2179;4061 2370;...
    3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
    3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
    2370 2975];                 %31個省會城市坐標
N=size(C,1);                    %TSP問題的規模,即城市數目,也就是基因數
D=zeros(N);                     %任意兩個城市距離間隔矩陣
%%%%%%%%%%%%%%%%%%%%%求任意兩個城市距離間隔矩陣%%%%%%%%%%%%%%%%%%%%%
for i=1:N
    for j=1:N
        D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;%計算距離
    end
end

NP=200;                          %種群規模
G=1000;                          %最大遺傳代數
f=zeros(NP,N);                   %用於存儲種群

F=[];                            %種群更新中間存儲
for i=1:NP
    f(i,:)=randperm(N);          %隨機生成初始種群
end
R=f(1,:);                        %存儲最優種群
len=zeros(NP,1);                 %存儲路徑長度
fitness=zeros(NP,1);             %存儲歸一化適應值


gen=0;
%%%%%%%%%%%%%%%%%%%%%%%%%遺傳算法循環%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while gen<G
    %%%%%%%%%%%%%%%%%%%%%計算路徑長度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:NP
        len(i,1)=D(f(i,N),f(i,1));	%終點到起點的距離
        for j=1:(N-1)
            len(i,1)=len(i,1)+D(f(i,j),f(i,j+1));	%路徑長度累加
        end
    end
    maxlen=max(len);              %最長路徑
    minlen=min(len);              %最短路徑
    %%%%%%%%%%%%%%%%%%%%%%%%%更新最短路徑%%%%%%%%%%%%%%%%%%%%%%%%%%
    rr=find(len==minlen);	%最短路徑位置數組
    R=f(rr(1,1),:);			%最短路徑
    %%%%%%%%%%%%%%%%%%%%%計算歸一化適應值%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:length(len)
        fitness(i,1)=(1-((len(i,1)-minlen)/(maxlen-minlen+0.001)));
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%選擇操作%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nn=0;
    for i=1:NP
        if fitness(i,1)>=rand  %適應度越高越有機會被選中
            nn=nn+1;         %被選中個體數加1
            F(nn,:)=f(i,:);  %用F記錄被選中的個體
        end
    end
    [aa,bb]=size(F);
    while aa<NP
        nnper=randperm(nn); %對選擇的群體里隨機挑2個個體
        A=F(nnper(1),:);    %隨機序列第一位對應的個體
        B=F(nnper(2),:);    %隨機序列第二位對應的個體
        %%%%%%%%%%%%%%%%%%%%%%%交叉操作%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        W=ceil(N/10);              	%交叉點個數31/10 = 4個
        p=unidrnd(N-W+1);          	%隨機選擇交叉范圍,從p到p+W
        for i=1:W    				%4個交叉點
        	%先保存交叉點所在值在對方的位置,然后交換交叉點值,再交換失去的值在對方的位置,即保證31個城市不會因為交叉導致沒走完
            x=find(A==B(p+i-1)); 	%B中交叉點所在值在A的位置
            y=find(B==A(p+i-1));
            temp=A(p+i-1);
            A(p+i-1)=B(p+i-1); 
            B(p+i-1)=temp;
            temp=A(x); 
            A(x)=B(y); 
            B(y)=temp;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%變異操作%%%%%%%%%%%%%%%%%%%%%%%%%
        p1=floor(1+N*rand());
        p2=floor(1+N*rand());
        while p1==p2				%生成兩個不同城市位置序號
            p1=floor(1+N*rand());
            p2=floor(1+N*rand());
        end
        tmp=A(p1); 					%在A與B中進行該位置上城市的交換,即發生變異
        A(p1)=A(p2); 
        A(p2)=tmp;
        tmp=B(p1); 
        B(p1)=B(p2); 
        B(p2)=tmp;
        F=[F;A;B];					%合並取出來並經過交叉,變異的個體到原選擇數組上
        [aa,bb]=size(F);
    end
    if aa>NP					 	
        F=F(1:NP,:);             	%保持種群規模為NP
    end
    f=F;                         	%更新種群
    f(1,:)=R;                    	%保留每代最優個體
    clear F;
    gen=gen+1
    Rlength(gen)=minlen;
    
end

figure
for i=1:N-1
    plot([C(R(i),1),C(R(i+1),1)],[C(R(i),2),C(R(i+1),2)],'bo-');
    hold on;
end
plot([C(R(N),1),C(R(1),1)],[C(R(N),2),C(R(1),2)],'ro-');
title(['優化最短距離:',num2str(minlen)]);
figure
plot(Rlength)
xlabel('迭代次數')
ylabel('目標函數值')
title('適應度進化曲線')

差分進化算法DE

差分演化算法是一種基於群體差異的啟發式隨機搜索算法,將問題的求解表示成"染色體"的適者生存過程,通過"染色體"群的一代代不斷進化,包括復制、交叉和變異等操作,最終收斂到"最適應環境"的個體,從而求得問題的最優解或滿意解。

關鍵參數

• 種群數量Np:一般來說,種群數量越大,種群的多樣性也就越大,尋優能力也就越強,但也會增加計算的難度。為了算法具有足夠的不同的變異向量,Np >= 4。一般取5D~10D。(D為變量維數)

• 變異算子F:變異算子 F∈[0, 2] 決定偏差向量的縮放比例。F過小,可能造成算法“早熟”,容易陷入局部最優;F過大,算法很難快速收斂到最優值。F = 0.5通常是一個較好的初始選擇。如果種群過早收斂,那么 F 或 Np 應該增大。

• 交叉算子CR:交叉算子 CR∈[0, 1]控制着一個試驗向量參數來自隨機選擇的變異向量而不是原來的向量的概率。CR越大,發生交叉的概率越大。CR的一個較好選擇是 0.1,但較大的 CR通常會加速收斂。為了嘗試獲得一個快速解,可以先嘗試 CR = 0.9 或 CR = 0.1。

主要流程

變異操作

DE/rand/1/bin

描述:

  1. 對需要變異的染色體 m,通過隨機值選擇 3 條互不相同且與染色體 m 也不相同的其他染色體 r1, r2, r3
  2. 利用變異算子對 r2, r3 的差異進行縮小,
  3. r2, r3與 r1 進行合並后替換到原染色體 m 上,完成變異過程。

DE/rand/1/bin + 自適應變異

描述:

在 DE/rand/1/bin 的基礎上,針對變異算子進行自適應計算。此自適應計算的方法可以基於變異代數,對於變異代數越高的變異操作,可以考慮減小變異算子,從而做到更好的收斂效果。

自適應算子計算

將變異算子中隨機選擇的三個個體Xb,Xm,Xw進行計算得變異種群數組

交叉

DE 交叉

描述:

同時采用隨機選擇和概率遺傳兩種方法,隨機選擇通過隨機選擇一位基因進行變異替換,概率遺傳通過概率值決定是否采用變異后的基因進行交叉。這樣可以確保子代至少有一個基因來自變異。

二項式交叉

選擇

貪婪准則:對原種群和變異后種群每個個體一對一進行判斷,取適應度更高的作為下一代種群

代碼

DE/rand/1/bin + 自適應變異

%%%%%%%%%%%%%%%%%差分進化算法求函數極值%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                            %清除所有變量
close all;                            %清圖
clc;                                  %清屏
NP=50;                                %個體數目
D=10;                                 %變量的維數
G=200;                                %最大進化代數
F0=0.4;                               %初始變異算子
CR=0.1;                               %交叉算子
Xs=20;                                %上限
Xx=-20;                               %下限
yz=10^-6;                             %閾值
MAXRUN = 10;        %獨立測試次數
%%%%%%%%%%%%%%%%%%%%%%%%%賦初值%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(D,NP);                        %初始種群
v=zeros(D,NP);                        %變異種群
u=zeros(D,NP);                        %選擇種群

for run = 1:MAXRUN
    %%%%%%%%%%%%%%%%一次獨立測試%%%%%%%%%%%%%%
    x=rand(D,NP)*(Xs-Xx)+Xx;       	%賦初值
   	%%%%%%%%%%%%%%%%%%%%計算目標函數%%%%%%%%%%%%%%%%%%%%
    for m=1:NP
        Ob(m)=func1(x(:,m));		%計算初始種群每位上的目標函數值
    end
    trace(1)=min(Ob);				%獲取最優值

    %%%%%%%%%%%%%%%%%%%%%%%差分進化循環%%%%%%%%%%%%%%%%%%%%%
    for gen=1:G
        %%%%%%%%%%%%%%%%%%%%%%變異操作%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%自適應變異算子%%%%%%%%%%%%%%%%%%%
        lamda=exp(1-G/(G+1-gen));
        F=F0*2^(lamda);				%自適應變異算子公式
        %%%%%%%%%%%%%%%%%r1,r2,r3和m互不相同%%%%%%%%%%%%%%%%
        for m=1:NP
            r1=randi([1,NP],1,1);
            while (r1==m)						%防止r1==m
                r1=randi([1,NP],1,1);
            end
            r2=randi([1,NP],1,1);
            while (r2==m)|(r2==r1)				%防止r2==r1,r2==m
                r2=randi([1,NP],1,1);
            end
            r3=randi([1,NP],1,1);
            while (r3==m)|(r3==r1)|(r3==r2)		%防止r3==r1,r3==r2,r3==m
                r3=randi([1,NP],1,1);
            end
            v(:,m)=x(:,r1)+F*(x(:,r2)-x(:,r3));	%變異向量公式3.3,通過對r1,r2,r3列操作后放入變異種群數組的第m列
        end
        %%%%%%%%%%%%%%%%%%%%%%交叉操作%%%%%%%%%%%%%%%%%%%%%%%
        r=randi([1,D],1,1);			%整個交叉操作的隨機值,確保交叉操作肯定會選擇變異種群其中一個作為新種群,隨機選擇法
        for n=1:D
            cr=rand(1);				%每個變量上的隨機值,概率遺傳法
            if (cr<=CR)|(n==r)		%u為新選擇種群,v為變異種群,x為初始種群
            %如果每個變量上的隨機值小於變異算子,或n==整個操作的隨機值,則進行交叉,即用當前序號變異種群作為選擇種群
                u(n,:)=v(n,:);
            else					%否則,則不進行交叉,即用當前序號初始種群作為選擇種群
                u(n,:)=x(n,:);
            end
        end
        %%%%%%%%%%%%%%%%%%%邊界條件的處理%%%%%%%%%%%%%%%%%%%%%
        for n=1:D
            for m=1:NP
                if (u(n,m)<Xx)|(u(n,m)>Xs)		%超過邊界的重新生成值
                    u(n,m)=rand*(Xs-Xx)+Xx;
                end
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%選擇操作%%%%%%%%%%%%%%%%%%%%%%%
        for m=1:NP
            Ob1(m)=func1(u(:,m));	%重新計算新選擇種群每位上的目標函數值
        end
        for m=1:NP
            if Ob1(m)<Ob(m)			%若小於之前的初始種群的該位上的目標函數值(更優),初始種群上該位被選擇種群替換
                x(:,m)=u(:,m);		%即是適者生存
            end
        end  
        for m=1:NP
            Ob(m)=func1(x(:,m));	%重新計算現在的初始種群每位上的目標函數值
        end
        trace(gen+1)=min(Ob);		%取最優值放入軌跡
        if min(Ob(m))<yz			%如果最優值小於閾值,則跳出
            break
        end
end
[SortOb,Index]=sort(Ob);
x=x(:,Index);
X=x(:,1);                          	%最優變量              
Y=min(Ob)                           %最優值  
end
%%%%%%%%%%%%%%%%%%%%%%%%%畫圖%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(trace);
xlabel('迭代次數')
ylabel('目標函數值')
title('適應度進化曲線')

免疫算法IA

免疫算法是模仿生物免疫機制,結合基因的進化機理,人工構造出的一種新型智能優化算法。采用群體搜索策略,通過迭代計算,最終以較大的概率得到問題的最優解。相比於其他算法,免疫算法利用自身產生多樣性和維持機制的特點,保證了群體多樣性,克服了‘早熟’問題,可以得到全局最優解。具有自適應性,隨機性並行性,全局收斂性,種群多樣性等特點。

生物免疫系統 免疫算法
抗原 優化問題
抗體(B細胞) 優化問題的可行解
親和度 可行解的質量
細胞活化 免疫選擇
細胞分化 個體克隆
親和度成熟 變異
克隆抑制 克隆抑制
動態維持平衡 種群刷新

關鍵參數

  • 抗體種群大小Np:保留了免疫細胞的多樣性,種群越大,算法全局搜索能力越好,但算法每代的計算量也就相應增大。一般取10~100,且200以內。
  • 免疫選擇比例:免疫選擇的抗體數量越多,產生抗體越多,搜索能力越強,但增加每代計算量。一般取0.1Np
    0.5Np
  • 抗體克隆擴增的倍數:決定了克隆擴增的細胞的數量,從而決定算法的搜索能力。數值越大,局部搜索能力越強,全局搜索能力也有一定提高,但是計算量增大。一般取5~10。
  • 種群刷新比例:細胞的淘汰和更新是產生抗體多樣性的重要機制,從而對免疫算法的全局搜索能力產生重要影響。一般不超過0.5Np。
  • 算子
    • 親和度評價算子:相當於遺傳算法中的適應度,與具體問題相關
    • 抗體濃度評價算子:表征抗體種群的多樣性好壞,濃度過高意味非常類似個體大量存在;不利於全局優化
    • 激勵度算子:對抗體質量的最終評價結果,通常親和度大,濃度低的抗體會得到較大激勵度
    • 克隆抑制算子:用於對經過變異后的克隆體進行再選擇,抑制親和度低的抗體(重新生成隨機新生種群-種群刷新),保留親和度高的抗體作為免疫種群,最后免疫種群與新生種群進行進行合並

主要流程

免疫選擇

前 NP / 2 個激勵度高的個體指針對 NP/2 個激勵度高的個體進行交叉和變異,目的是試圖產生激勵度更高的個體。隨機生成指針對另一半的個體產生采用隨機生成的方式,保證了種群中不斷有新個體的加入。

抗體間親和度計算

對於實數編碼,親和度通常可以使用抗體向量間的歐氏距離來計算

抗體濃度計算

抗體濃度公式

其中N為種群規模,S(abi, abj)表示抗體間的相似度,具體表示如下

相似度檢查

若小於相似度閾值,則設為1

激勵度計算

激勵度算子就是抗體的綜合評價,計算方式如下,二者取其一,其本質目的是為了共同考慮親和度和濃度,從而篩選下一代的抗體

變異操作

變異源鄰域

描述:

在需要變異的個體染色體中,通過添加一個擾動,從而使其偏離原來的位置落入原個體相鄰的另外一個位置中,從而實現變異源領域的搜索。

實數編碼算法變異計算

交叉操作

克隆抑制

描述:

對於克隆產生的所有個體進行親和度的計算,然后保留親和度最高的個體,從而實現交叉操作。

例子:

  1. 原個體親和度為 10
  2. 克隆后新個體親和度分別為 [12,6,34,43,2,15,3,8,19,22]
  3. 將10個新個體和原個體合並后取親和度最高的個體,即親和度為 43 的個體

代碼

%%%%%%%%%%%%%%%%%免疫算法求函數極值%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                                %清除所有變量
close all;                                %清圖
clc;                                      %清屏
D=10;                                     %免疫個體維數
NP=100;                                   %免疫個體數目
Xs=20;                                    %取值上限
Xx=-20;                                   %取值下限
G=500;                                    %最大免疫代數
pm=0.7;                                   %變異概率
alfa=1;                                   %激勵度系數
belta=1;                                  %激勵度系數   
detas=0.2;                                %相似度閾值
gen=0;                                    %免疫代數
Ncl=10;                                   %克隆個數
deta0=1*Xs;                               %鄰域范圍初值
%%%%%%%%%%%%%%%%%%%%%%%初始種群%%%%%%%%%%%%%%%%%%%%%%%%
f=rand(D,NP)*(Xs-Xx)+Xx;	%隨機生成初始值數組10行100列
for np=1:NP	
    FIT(np)=func1(f(:,np));
end
%%%%%%%%%%%%%%%%%計算個體濃度和激勵度%%%%%%%%%%%%%%%%%%%
for np=1:NP
    for j=1:NP     
        nd(j)=sqrt(sum((f(:,np)-f(:,j)).^2));%基於歐式距離的抗體間親和度計算
        %nd一開始存儲抗體間親和度,后來計算出抗體間相似度,從而取平均值獲得抗體濃度,最后通過
        if nd(j)<detas	%如果小於相似度閾值,則說明相似,設為1
            nd(j)=1;
        else
            nd(j)=0;	%如果大於等於相似度閾值,則說明不相似,設為0
        end
    end
    ND(np)=sum(nd)/NP;	%取平均值得抗體濃度
end
FIT =  alfa*FIT- belta*ND;	%4.6 抗體abi激勵度=激勵度系數ax初始激勵度-激勵度系數bx抗體濃度
%%%%%%%%%%%%%%%%%%%激勵度按升序排列%%%%%%%%%%%%%%%%%%%%%%
[SortFIT,Index]=sort(FIT);
Sortf=f(:,Index);
%%%%%%%%%%%%%%%%%%%%%%%%免疫循環%%%%%%%%%%%%%%%%%%%%%%%%
while gen<G
    for i=1:NP/2
        %%%%%%%%選激勵度前NP/2個體進行免疫操作%%%%%%%%%%%
        a=Sortf(:,i);		%取第i列
        Na=repmat(a,1,Ncl);	%生成重復第i列的1xNcl的矩陣
        deta=deta0/gen;		%鄰域范圍初值/當前循環代數=擾動因子,代數越大,擾動越低,減小誤差
        for j=1:Ncl			%迭代克隆個體數
            for ii=1:D
                %%%%%%%%%%%%%%%%%變異%%%%%%%%%%%%%%%%%%%
                if rand<pm	%如果小於變異率,則發生變異
                    Na(ii,j)=Na(ii,j)+(rand-0.5)*deta;	%變異算子
                end
                %%%%%%%%%%%%%%邊界條件處理%%%%%%%%%%%%%%%
                if (Na(ii,j)>Xs)  |  (Na(ii,j)<Xx)
                    Na(ii,j)=rand * (Xs-Xx)+Xx;		%如果超出邊界則重新生成
                end
            end
        end
        Na(:,1)=Sortf(:,i);             %保留克隆源個體
        %%%%%%%%%%克隆抑制,保留親和度最高的個體%%%%%%%%%%
        for j=1:Ncl
            NaFIT(j)=func1(Na(:,j));	%親和度計算,得數組
        end
        [NaSortFIT,Index]=sort(NaFIT);	%親和度排序
        aFIT(i)=NaSortFIT(1);			%取函數最小值,即親和度最高值
        NaSortf=Na(:,Index);			
        af(:,i)=NaSortf(:,1);			%親和度最高的個體
    end 
    %%%%%%%%%%%%%%%%%%%%免疫種群激勵度%%%%%%%%%%%%%%%%%%%
    for np=1:NP/2
        for j=1:NP/2
            nda(j)=sqrt(sum((af(:,np)-af(:,j)).^2));         
            if nda(j)<detas
                nda(j)=1;
            else
                nda(j)=0;
            end
        end
        aND(np)=sum(nda)/NP/2;
    end
    aFIT =  alfa*aFIT-  belta*aND;
    %%%%%%%%%%%%%%%%%%%%%%%種群刷新%%%%%%%%%%%%%%%%%%%%%%%
    bf=rand(D,NP/2)*(Xs-Xx)+Xx;		%生成D行NP/2列在取值范圍內的新種群
    for np=1:NP/2
        bFIT(np)=func1(bf(:,np));	%計算新種群函數評價
    end
    %%%%%%%%%%%%%%%%%%%新生成種群激勵度%%%%%%%%%%%%%%%%%%%%
    for np=1:NP/2
        for j=1:NP/2
            ndc(j)=sqrt(sum((bf(:,np)-bf(:,j)).^2));
            if ndc(j)<detas
                ndc(j)=1;
            else
                ndc(j)=0;
            end
        end
        bND(np)=sum(ndc)/NP/2;
    end
    bFIT =  alfa*bFIT-  belta*bND;
    %%%%%%%%%%%%%%免疫種群與新生種群合並%%%%%%%%%%%%%%%%%%%
    f1=[af,bf];					%新舊種群合並
    FIT1=[aFIT,bFIT];			%新舊種群激勵度數組合並
    [SortFIT,Index]=sort(FIT1);	%新舊種群激勵度合並數組排序
    Sortf=f1(:,Index);			%取激勵度最高的一列,即這一代最優解
    gen=gen+1;					%免疫代數+1
    trace(gen)=func1(Sortf(:,1));	%加入最優解軌跡數組
end
%%%%%%%%%%%%%%%%%%%%%%%輸出優化結果%%%%%%%%%%%%%%%%%%%%%%%%
Bestf=Sortf(:,1);                 %最優變量
trace(end);                       %最優值
figure,plot(trace)
xlabel('迭代次數')
ylabel('目標函數值')
title('親和度進化曲線')

蟻群算法ACO

蟻群算法是一種仿生學算法,是由自然界中螞蟻覓食的行為而啟發的。

在自然界中,螞蟻覓食過程中,蟻群總能夠按照尋找到一條從蟻巢和食物源的最優路徑。螞蟻在尋找食物的過程中往往是隨機選擇路徑的,但它們能感知當前地面上的信息素濃度, 並傾向於往信息素濃度高的方向行進。信息素由螞蟻自身釋放,是實現蟻群內間接通信的物質。

由於較短路徑上螞蟻的往返時間比較短,單位時間內經過該路徑的螞蟻多,所以信息素的積累速度比較長路徑快。因此,當后續螞蟻在路口時,就能感知先前螞蟻留下的信息,並 傾向於選擇一條較短的路徑前行。

這種正反饋機制使得越來越多的螞蟻在巢穴與食物之間的 最短路徑上行進。由於其他路徑上的信息素會隨着時間蒸發,最終所有的螞蟻都在最優路徑上行進。

  • 構建解
    • 狀態轉移概率
      • 啟發式信息
        • 啟發因子:城市間距離的倒數
      • 信息素
  • 信息素
    • 蒸發
    • 增量
      • 最優解
      • 走過部分解

關鍵參數

  • 信息素啟發式因子α:代表信息量對是否選擇當前路徑的影響程度,反應螞蟻在運動過程中所積累的信息素在知道蟻群搜索中的相對重要程度。一般取α ∈[1, 4]。
  • 期望啟發因子β:表示在搜索時路徑上的信息素在指導螞蟻選擇路徑時的向導性,反映蟻群在搜索最優路徑的過程中的先驗性和確定性因素的作用強度。一般取β ∈[3, 5]。β越大,螞蟻在某個局部點上選擇局部最短路徑可能性越大,算法收斂速度越快,但隨機性越弱,容易陷入局部最優解
  • 信息蒸發系數p:p ∈[0,1],表示信息素的蒸發程度,反映了螞蟻群體中個體之間相互影響的強弱。p過小時,影響算法的隨機性能和全局搜索能力;p過大是,降低算法的收斂速度。
  • 信息素強度Q:代表計算信息素增量的一個取值
  • 螞蟻數目m:螞蟻數量增多,提高算法的全局搜索能力和穩定性,但正反饋作用不明顯,收斂速度減慢;螞蟻數量減少,收斂速度加快,但算法全局性能降低,穩定性差,容易出現過早停滯現象。一般取10~50。

主要流程

選擇

采用的是輪盤賭的方式,其概率的計算是基於信息素和距離權重的,公式如下(![img](file:///C:\Users\79304\AppData\Local\Temp\ksohtml1764\wps2.png) 為信息素,![img](file:///C:\Users\79304\AppData\Local\Temp\ksohtml1764\wps3.png) 為距離權重)

信息素更新

信息素增量計算公式,其中 Q 信息素強度系數,![img](file:///C:\Users\79304\AppData\Local\Temp\ksohtml1764\wps4.png) 為第 k 只螞蟻在本輪周游中所走過的路徑的長度

信息素蒸發量計算公式,其中 ![img](file:///C:\Users\79304\AppData\Local\Temp\ksohtml1764\wps5.png) 表示路徑上信息素的蒸發系數,![img](file:///C:\Users\79304\AppData\Local\Temp\ksohtml1764\wps6.png)表示信息素的持久性系數

代碼

蟻群算法解決旅行商問題

%%%%%%%%%%%%%%%%%%%%蟻群算法解決TSP問題%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                %清除所有變量
close all;                %清圖
clc;                      %清屏
m=50;                     %螞蟻個數
Alpha=1;                  %信息素重要程度參數              
Beta=5;                   %啟發式因子重要程度參數
Rho=0.1;                  %信息素蒸發系數
G_max=200;                %最大迭代次數
Q=100;                    %信息素增加強度系數
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
    3238 1229;4196 1044;4312  790;4386  570;3007 1970;2562 1756;...
    2788 1491;2381 1676;1332  695;3715 1678;3918 2179;4061 2370;...
    3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
    3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
    2370 2975];                 %31個省會城市坐標
%%%%%%%%%%%%%%%%%%%%%%%%第一步:變量初始化%%%%%%%%%%%%%%%%%%%%%%%%
n=size(C,1);              	%n表示問題的規模(城市個數)
D=zeros(n,n);             	%D表示兩個城市距離間隔矩陣
for i=1:n
    for j=1:n
        if i~=j				%如果不是同一個城市
            D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;	%計算兩個城市距離
        else
            D(i,j)=eps;		%設為一個數可以分辨的最小精度
        end
        D(j,i)=D(i,j);		%更新兩個城市距離
    end
end
Eta=1./D;                    %Eta為啟發因子,這里設為距離的倒數
Tau=ones(n,n);               %Tau為信息素矩陣
Tabu=zeros(m,n);             %存儲並記錄路徑的生成
NC=1;                        %迭代計數器
R_best=zeros(G_max,n);       %各代最佳路線
L_best=inf.*ones(G_max,1);   %各代最佳路線的長度
figure(1);%優化解
while NC<=G_max            
    %%%%%%%%%%%%%%%%%%第二步:將m只螞蟻放到n個城市上%%%%%%%%%%%%%%%%
    Randpos=[];
    for i=1:(ceil(m/n))					%求得每個城市放m/n只螞蟻,總共需要幾次迭代
        Randpos=[Randpos,randperm(n)];	%每次迭代取n只螞蟻隨機去往n個城市
    end
    Tabu(:,1)=(Randpos(1,1:m))'; 		%最后取存儲m只螞蟻去往的城市序號的數組
    %%%%%第三步:m只螞蟻按概率函數選擇下一座城市,完成各自的周游%%%%%%
    for j=2:n
        for i=1:m
            visited=Tabu(i,1:(j-1));  	%已訪問的城市數組
            J=zeros(1,(n-j+1));       	%待訪問的城市數組,未初始化,全為0
            P=J;                      	%待訪問城市的選擇概率分布
            Jc=1;						%待訪問的城市數組的遍歷變量
            for k=1:n
                if length(find(visited==k))==0	%find查找已訪問城市為k的,生成序號數組,並計算有幾個(length)
                    J(Jc)=k;					%若為0,則待訪問,將城市序號用來初始化待訪問的城市數組
                    Jc=Jc+1;					%遍歷變量++
                end
            end
            %%%%%%%%%%%%%%%%%%計算待選城市的概率分布%%%%%%%%%%%%%%%%
            for k=1:length(J)
                P(k)=(Tau(visited(end),J(k))^Alpha)...
                    *(Eta(visited(end),J(k))^Beta);
            end
            P=P/(sum(P));	%狀態轉移概率公式計算5.1
            %%%%%%%%%%%%%%%%按概率原則選取下一個城市%%%%%%%%%%%%%%%%
            Pcum=cumsum(P);				%輪盤賭概率累加
            Select=find(Pcum>=rand);	%隨機生成一個概率,查找概率累加數組里大於等於該隨機概率的數,生成數組Select
            to_visit=J(Select(1));		%取Select數組第一個數作為下一個城市
            Tabu(i,j)=to_visit;
        end
    end
    if NC>=2
        Tabu(1,:)=R_best(NC-1,:);
    end
    %%%%%%%%%%%%%%%%%%%第四步:記錄本次迭代最佳路線%%%%%%%%%%%%%%%%%%
    L=zeros(m,1);
    for i=1:m
        R=Tabu(i,:);					%第i行路線放入R
        for j=1:(n-1)					
            L(i)=L(i)+D(R(j),R(j+1));	%遍歷n個城市計算第i行路線長度
        end
        L(i)=L(i)+D(R(1),R(n));			%最后補上終點到起點的距離
    end
    L_best(NC)=min(L);					%這一代最佳路徑的長度記錄
    pos=find(L==L_best(NC));			%查找哪一行(路線)長度==最佳路徑長度
    R_best(NC,:)=Tabu(pos(1),:);		%這一代最佳路徑的記錄
    %%%%%%%%%%%%%%%%%%%%%%%%%第五步:更新信息素%%%%%%%%%%%%%%%%%%%%%%
    Delta_Tau=zeros(n,n);
    for i=1:m
        for j=1:(n-1)
            Delta_Tau(Tabu(i,j),Tabu(i,j+1))=...		%ant-cycle模型
                Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);%Q/L(i)=信息素增加強度系數/本次周游所走路徑長度=迭代新增信息素
        end
        Delta_Tau(Tabu(i,n),Tabu(i,1))=...
            Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);
    end
    Tau=(1-Rho).*Tau+Delta_Tau;		%5.2信息素公式=保留下來的信息素+迭代新增信息素,Rho代表蒸發系數,1-Rho代表持久系數
    %%%%%%%%%%%%%%%%%%%%%%%第六步:禁忌表清零%%%%%%%%%%%%%%%%%%%%%%
    Tabu=zeros(m,n);
    %%%%%%%%%%%%%%%%%%%%%%%%%歷代最優路線%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:n-1
        plot([ C(R_best(NC,i),1), C(R_best(NC,i+1),1)],...
            [C(R_best(NC,i),2), C(R_best(NC,i+1),2)],'bo-');	%畫路線
        hold on;
    end
    plot([C(R_best(NC,n),1), C(R_best(NC,1),1)],...
        [C(R_best(NC,n),2), C(R_best(NC,1),2)],'ro-');  
    title(['優化最短距離:',num2str(L_best(NC))]);
    hold off;
    pause(0.005);
    NC=NC+1;    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%第七步:輸出結果%%%%%%%%%%%%%%%%%%%%%%%%%%
Pos=find(L_best==min(L_best));
Shortest_Route=R_best(Pos(1),:);            %最佳路線
Shortest_Length=L_best(Pos(1));             %最佳路線長度
figure(2),
plot(L_best)
xlabel('迭代次數')
ylabel('目標函數值')
title('適應度進化曲線')

粒子群算法PSO

粒子群算法模擬鳥群的捕食行為,每個優化問題的解都是搜索空間中的一只鳥。我們稱之為“粒子”。所有的粒子都有一個由被優化的函數決定的適應值(fitnessvalue),每個粒子還有一個速度決定他們飛翔的方向和距離。然后粒子們就追隨當前的最優粒子在解空間中搜索。

關鍵參數

  • 粒子種群規模N:粒子數目越大,算法搜索的空間范圍就越大,也就更容易發現全局最優解,但算法運行時間也就越長。一般10個粒子已經可以取得比較好的結構。對於比較難的問題或特定的問題,粒子的數量可以取到100或200。
  • 慣性權重w:代表對粒子當前速度繼承的多少,用來控制算法的開發和探索能力。權重較大時,全局尋優能力較強,局部尋優能力較弱;權重較小時,全局尋優能力較弱,局部尋優能力較強。一般取0.8~1.2。
  • 加速常數c1和c2:分別調節向 pbest 和 gbest方向飛行的最大步長,他們分別決定粒子個體經驗和群體經驗對粒子運行軌跡的影響。通常可以取c1 = c2 = 1.5。

主要流程

PSO初始化為一群隨機粒子(隨機解),然后通過迭代找到最優解,在每一次迭代中,粒子通過跟蹤兩個“極值”來更新自己。

第一個就是粒子本身所找到的最優解,這個解叫做個體極值pBest,

另一個極值是整個種群找到的最優解,這個極值是全局極值gBest。另外也可以不用整個種群而只是用其中一部分最優粒子的鄰居,那么在所有鄰居中的極值就是局部極值。

  • 粒子的速度及位置更新的方式如下:

[公式]

  1. 其中 [公式]
  2. 粒子的速度更新公式由三部分組成:粒子先前速度+個體最優+全局最優
  3. [公式] 是一個非負數,稱為慣性因子,對算法的收斂起到很大的作用,其值越大,粒子飛躍的范圍就越廣,更容易找到全局最優,但是也會錯失局部搜尋的能力。
  4. [公式] 分別為局部和全局最優位置。
  5. 加速常數 [公式] 也是非負常數,是調整局部最優值和全局最優值權重的參數,如果前者為0說明搜尋過程中沒有自身經驗只有社會經驗,容易陷入局部最優解;若后者為0,即只有社會經驗,沒有自身經驗,常常會陷入局部最優解中,不能飛越該局部最優區域。
  6. [公式] 是[0,1]范圍之內的隨機數, [公式] 是約束因子,目的是控制速度的權重。

代碼

%%%%%%%%%%%%%%%%%粒子群算法求函數極值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                  %清除所有變量
close all;                  %清圖
clc;                        %清屏
global fes;                 %fes統計函數計算次數
N = 100;                    %群體粒子個數
D = 10;                     %粒子維數
c1 = 1.5;                   %加速常數1
c2 = 1.5;                   %加速常數2
w = 0.8;                    %慣性權重
a = 1.0;					%約束因子
Xmax = 20;                  %位置最大值
Xmin = -20;                 %位置最小值
Vmax = 10;                  %速度最大值
Vmin = -10;                 %速度最小值

MAXRUN = 20;                %獨立測試次數
MAXFES = 10^5;              %最大函數迭代次數
error = 10^-6;              %誤差閾值
success = 0;

for run = 1:MAXRUN
    fes = 0;   
    %%%%%%%%%%%%%%%%初始化種群個體(限定位置和速度)%%%%%%%%%%%%%%%%
    x = rand(N, D) * (Xmax - Xmin) + Xmin;
    v = rand(N, D) * (Vmax - Vmin) + Vmin;
    %%%%%%%%%%%%%%%%%%初始化個體最優位置和最優值%%%%%%%%%%%%%%%%%%%
    p = x;
    pbest = ones(N,1);
    for i = 1:N
        pbest(i) = func(x(i, :));
    end
    %%%%%%%%%%%%%%%%%%%初始化全局最優位置和最優值%%%%%%%%%%%%%%%%%%
    g = ones(1, D);
    gbest = inf;				%初始設為正無窮
    for i = 1:N
        if(pbest(i) < gbest)	%如果個體最優小於全局最優
            g = p(i, :);		%更新全局最優數組
            gbest = pbest(i);	%更新全局最優值
        end
    end
    
    gen = 1;
    %%%%%%%%%%%按照公式依次迭代直到滿足精度或者迭代次數%%%%%%%%%%%%%
    while fes <= MAXFES
        for j = 1:N
            %%%%%%%%%%%%%%更新個體最優位置和最優值%%%%%%%%%%%%%%%%%
            temp = func(x(j, :));
            if (temp < pbest(j))
                p(j, :) = x(j, :);
                pbest(j) = temp;
            end
            %%%%%%%%%%%%%%%%更新全局最優位置和最優值%%%%%%%%%%%%%%%
            if(pbest(j) < gbest)
                g = p(j, :);
                gbest = pbest(j);
            end
            %%%%%%%%%%%%%%%%%更新位置和速度值%%%%%%%%%%%%%%%%%%%%%
            v(j, :) = w * v(j, :) + c1 * rand * (p(j, :) - x(j, :))...	%粒子的速度更新公式
                + c2 * rand * (g - x(j, :));							
            x(j, :) = x(j, :) + a*v(j, :);								%粒子的位置更新公式
            %%%%%%%%%%%%%%%%%%%%邊界條件處理%%%%%%%%%%%%%%%%%%%%%%
            for ii = 1:D
                if (v(j, ii) > Vmax) | (v(j, ii) < Vmin)		%對粒子的速度,位置進行邊界處理
                    v(j, ii) = rand * (Vmax - Vmin) + Vmin;
                end
                if (x(j, ii) > Xmax) | (x(j,ii) < Xmin)
                    x(j, ii) = rand * (Xmax - Xmin) + Xmin;
                end
            end
        end
        %%%%%%%%%%%%%%%%%%%%記錄歷代全局最優值%%%%%%%%%%%%%%%%%%%%%
        trace(gen) = gbest;
        gen = gen + 1;
    end 
end

%%%%%%%%%%%%%%%%%%%%%%%%%畫圖%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(trace)
xlabel('迭代次數');
ylabel('適應度值');
title('適應度進化曲線')

模擬退火算法

模擬退火算法來源於固體退火原理,將固體加溫至充分高,再讓其徐徐冷卻,加溫時,固體內部粒子隨溫升變為無序狀,內能增大,而徐徐冷卻時粒子漸趨有序,在每個溫度都達到平衡態,最后在常溫時達到基態,內能減為最小,即所求值

能量即目標函數,能量最低態即最優解,

溫度是Metropolis算法的一個重要控制參數,開始時T大,可以接受較差惡化解;隨着T減小,只能接受較好的惡化解了

根據Metropolis准則,粒子在溫度T時趨於平衡的概率為 exp(-ΔE/(kT)) ,其中 E 為溫度 T 時的內能,ΔE為其改變數, kBoltzmann 常數。 Metropolis 准則常表示為

主要流程

  1. 系統從一個能量狀態變化到另一個能量狀態時(根據當前溫度產生新解:鄰域隨機選取+邊界處理),相應的能量從Eold變化到Enew
  2. 若新狀態是全局最優,則更新全局最優解,保留上一個最優解
  3. Metropolis過程
    1. 若狀態向下(局部最優),則接受新解
    2. 若狀態向上(全局搜索),則以一定概率接受(exp(-ΔE/(kT))>rand)
  4. 在系統每次變化時,T也在冷卻,T(n+1)=KxT(n)。當T趨向於0(設置一個小值,如0.001)時,求得問題整體最優解

禁忌搜索算法

基本思想:

采用鄰域選優的搜索方法,為了逃離局部最優解,算法必須能夠接受劣解,也就是每一次得到的解不一定優於原來的解。

但是,一旦接受了劣解,算法迭代即可能陷入循環。為了避免循環,算法將最近接受的一些移動放在禁忌表中,在以后的迭代中加以禁止。即只有不再禁忌表中的較好解(可能比當前解差)才能接受作為下一代迭代的初始解。

隨着迭代的進行,禁忌表不斷更新,經過一定的迭代次數后,最早進入禁忌表的移動就從禁忌表中解禁退出。

主要流程

算法記憶匯總

  • 遺傳算法
    • 最早是由美國的 John Holland於20世紀70年代提出
    • 能有效求解NP問題(所有的非確定性多項式時間可解的判定問題構成NP類問題),以及非線性,多峰函數優化和多目標優化問題
  • 差分進化算法
    • Storn和Price於1995年首次提出
  • 免疫算法
    • 1958年澳大利亞學者Burnet率先提出克隆選擇原理
    • 1973年Jerne提出免疫系統的數學框架
  • 蟻群算法
    • Marco Dorigo於1992年在他的博士論文中提出
    • 蟻群算法是一種用來尋找優化路徑的概率型算法,靈感來源於螞蟻在尋找食物過程中發現路徑的行為
  • 粒子群算法
    • Eberhart博士和kennedy博士發明
  • 模擬退火算法
    • 最早的思想是由Metropolis等人於1953年提出。1983 年,Kirkpatrick 等成功地將退火思想引入到組合優化領域。
    • 模擬退火算法從某一較高初溫出發,伴隨溫度參數的不斷下降,結合概率突跳特性在解空間中隨機尋找目標函數目標函數)的全局最優解,即在局部最優解能概率性地跳出並最終趨於全局最優。
  • 禁忌搜索算法
    • Glover教授於1986年在一篇論文中首次提出
    • 從一個初始可行解出發,選擇一系列的特定搜索方向(移動)作為試探,選擇實現讓特定的目標函數值變化最多的移動。為了避免陷入局部最優解,TS搜索中采用了一種靈活的“記憶”技術,對已經進行的優化過程進行記錄和選擇,指導下一步的搜索方向,這就是Tabu表的建立。
  • 神經網絡算法
    • 最早由McCullochPitts提出形式神經元數學模型
    • Hopfield提出具有聯想記憶功能的Hopfield神經網絡,標志神經網絡取得重大進展
    • 用大量的簡單計算單元(神經元)構成非線性系統,在一定程度上模擬了大腦的信息處理、存儲和檢索等功能。
    • 常見網絡結構
      • 前饋神經網絡
      • 反饋神經網絡
        • BP網絡的誤差反向后傳學習算法,是最常用的神經網絡算法。它利用輸出后的誤差來估計輸出層的直接前導層誤差,再利用這個誤差更新前一層的誤差,如此反傳下去獲得所有其他各層的誤差估計。
      • 自組織網絡

測試輸出匯總

  • tall:運行計時

  • mean:平均

  • gen:代數

  • fes:函數評價次數

  • trace:最優值

注意輸出信息最終應該包括:

  • 測試的函數名Fname算法參數取值,如D, Pc, Pm, NP,MAXG, MAXRUN等

  • 表3最優函數值bestfun, 最差函數值worstfun,平均值meanfun,標准差stdfun,算法總體平均用時meantotaltimes(s)

  • 表4首次獲得至今最優解值對應的平均代數meanbestGen,對應的函數評價次數meanbestFEs,對應的用時meanbesttime(s),統計得到的成功次數success%

//單次測試F1run
//每次運行記錄F1totoal
//匯總測試total求最優值,平均值,平均用時


global fes                 %fes統計函數計算次數
MAXRUN = 10;        %獨立測試次數
funfile = fopen(‘F1total.txt’,‘a’);            %w 改為a表示添加(w表示覆蓋原有內容)


tallrecord = [];		      %用於記錄每次獨立測試的運行用時
bestfrecord = [];                             %用於記錄每次獨立測試得到的最優函數值
allfile = fopen(‘total.txt’,‘a’);            %w 改為a表示添加(w表示覆蓋原有內容)
for run = 1:MAXRUN
       %....忽略一段代碼…..  
       fprintf('第%d次運行\t第%d代\t%f\n',run,0,trace(1));
       onerunfile = fopen([‘F1_run’ num2str(run) ‘.txt’],‘w’);   %建立第run次運行的記錄文檔
       fprintf(onerunfile,‘%d\t%d\t%g\r\n’,0,fes,trace(1));       %寫入第0代種群結
       
       
       t0 = clock;	 %記錄當前時刻
       tall = 0;	 %統計計算總用時(排除寫文件用時)
       for gen=1:G       %差分進化循環
	   		%....忽略好多段代碼…..     
       		tall = tall + etime(clock,t0);      %累加計算用時(單位秒,排除寫文件用時)
            fprintf(onerunfile,'%d\t%d\t%g\r\n',gen,fes,trace(gen+1)); %寫入第gen代結果
            t0 = clock;      %記錄當前時刻

       end
       fclose(onerunfile); 
       tallrecord = [tallrecord,tall];     %添加第run次獨立測試的運行用時
       bestfrecord = [bestfrecord, trace(gen+1)]; %添加第run次獨立測試得到的最優函數值
       
       fprintf(funfile,‘%d\t%g\t%g\r\n’,run,trace(gen+1),tall); %輸出第幾次運行最優值和用時
End
fclose(funfile);
fprintf(allfile,'%s\t%g\t%g\t%g\r\n','F1',min(bestfrecord),mean(bestfrecord),mean(tallrecord));
%輸出匯總信息:【函數名,多次運行最優值,平均值,平均用時】


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM