粒子群優化演算法對BP神經網路優化 Matlab實現

未名w發表於2020-05-27

1、粒子群優化演算法

粒子群演算法(particle swarm optimization,PSO)由Kennedy和Eberhart在1995年提出,該演算法模擬鳥叢集飛行覓食的行為,鳥之間通過集體的協作使群體達到最優目的,是一種基於 Swarm Inteligence的優化方法。同遺傳演算法類似,也是一種基於群體疊代的,但並沒有遺傳演算法用的交叉以及變異,而是粒子在解空間追隨最優的粒子進行搜尋。PSO的優勢在於簡單容易實現同時又有深刻的智慧背景,既適合科學研究,又特別適合工程應用,並且沒有許多引數需要調整。

粒子群優化演算法源於1987年Reynolds對鳥群社會系統boids的模擬研究,boids是一個CAS。在boids中,一群鳥在空中飛行,每個鳥遵守以下三條規則:
1)避免與相鄰的鳥發生碰撞衝突;
2)儘量與自己周圍的鳥在速度上保持協調和一致;
3)儘量試圖向自己所認為的群體中靠近。
僅通過使用這三條規則,boids系統就出現非常逼真的群體聚集行為,鳥成群地在空中飛行,當遇到障礙時它們會分開繞行而過,隨後又會重新形成群體。
Reynolds僅僅將其作為CAS的一個例項作模擬研究,而並未將它用於優化計算中 。
Kennedy和Eberhart在中加入了一個特定點,定義為食物,鳥根據周圍鳥的覓食行為來尋找食物。他們的初衷是希望通過這種模型來模擬鳥群尋找食源的現象,然而實驗結果卻揭示這個模擬模型中蘊涵著很強的優化能力,尤其是在多維空間尋優中。
PSO中,每個優化問題的解都是搜尋空間中的一隻鳥。稱之為“粒子(Particle)”。所有的粒子都有一個由被優化的函式決定的適應值,每個粒子還有一個速度決定他們飛翔的方向和距離。然後粒子們就追隨當前的最優粒子在解空間中搜尋.
PSO 初始化為一群隨機粒子。然後通過疊代找到最優解。在每一次疊代中,粒子通過跟蹤兩個"極值"來更新自己。第一個就是粒子本身所找到的最優解。這個解叫做個體極值pBest. 另一個極值是整個種群目前找到的最優解。這個極值是全域性極值gBest。另外,也可以不用整個種群而只是用其中一部分的鄰居。
PSO演算法數學表示如下:
設搜尋空間為D維,總粒子數為n。第i個粒子位置表示為向量Xi=( xi1, xi2,…, xiD );第i個粒子 “飛行”歷史中的過去最優位置(即該位置對應解最優)為Pi=( pi1,pi2,…,piD ),其中第g個粒子的過去最優位置Pg為所有Pi ( i=1, …,n)中的最優;第i個粒子的位置變化率(速度)為向量Vi=(vi1, vi2,…, viD)。每個粒子的位置按如下公式進行變化(“飛行”):

 

 

其中,C1,C2為正常數,稱為加速因子;rand( )為[0,1]之間的隨機數;w稱慣性因子,w較大適於對解空間進行大範圍探查(exploration),w較小適於進行小範圍開挖(exploitation)。第d(1≤d≤D)維的位置變化範圍為[-XMAXd , XMAXd],速度變化範圍為[-VMAXd , VMAXd],迭代中若位置和速度超過邊界範圍則取邊界值。

粒子群初始位置和速度隨機產生,然後按公式(1)(2)進行迭代,直至找到滿意的解。目前,常用的粒子群演算法將全體粒子群(Global)分成若干個有部分粒子重疊的相鄰子群,每個粒子根據子群(Local)內歷史最優Pl調整位置,即公式(2)中Pgd換為Pld。

 

2、粒子群優化BP神經網路

神經網路與粒子群演算法的結合主要有兩種方式:一是利用粒子群演算法的全域性搜尋能力來優化神經網路的拓撲結構、連線權值和閾值,將粒子群演算法良好的全域性尋優能力與BP演算法良好的區域性尋優能力相結合,以提高神經網路的泛化能力和學習效能,從而改進神經網路的整體搜尋效率;二是將神經網路嵌入到粒子群演算法當中,利用神經網路良好的學習效能來改進粒子群演算法的優化效能,以提高粒子群演算法的收斂速度,減少計算的工作量。

3、粒子群演算法生成BP神經網路的權值和閾值,Matlab程式設計

                            
function [W1,W2,B1,B2]=pso(HiddenNum,InDim,OutDim)

% clc;clear 
% tic;                              %程式執行計時
E0=0.001;                        %允許誤差
MaxNum=100;                    %粒子最大迭代次數
narvs=InDim*HiddenNum+1+HiddenNum+HiddenNum;                         %目標函式的自變數個數=51
particlesize=30;                    %粒子群規模=200
c1=2;                            %每個粒子的個體學習因子,也稱為加速常數
c2=2;                            %每個粒子的社會學習因子,也稱為加速常數
w=0.6;                           %慣性因子
vmax=0.8;                        %粒子的最大飛翔速度
x=-5+10*rand(particlesize,narvs);     %粒子所在的位置
v=2*rand(particlesize,narvs);         %粒子的飛翔速度
%用inline定義適應度函式以便將子函式檔案與主程式檔案放在一起,
%目標函式是:y=1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))
%inline命令定義適應度函式如下:
fitness=inline('1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');
%inline定義的適應度函式會使程式執行速度大大降低
for i=1:particlesize%粒子群規模
    for j=1:narvs %目標函式自變數個數
        f(i)=fitness(x(i,j));
    end
end
personalbest_x=x;%30*3
personalbest_faval=f;%30*1 但是是
[globalbest_faval i]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
k=1;
while k<=MaxNum
    for i=1:particlesize
        for j=1:narvs
            f(i)=fitness(x(i,j));
        end
        if f(i)<personalbest_faval(i) %判斷當前位置是否是歷史上最佳位置
            personalbest_faval(i)=f(i);
            personalbest_x(i,:)=x(i,:);
        end
    end
    [globalbest_faval i]=min(personalbest_faval);
    globalbest_x=personalbest_x(i,:);
    for i=1:particlesize %更新粒子群裡每個個體的最新位置
        v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...
            +c2*rand*(globalbest_x-x(i,:));
        for j=1:narvs    %判斷粒子的飛翔速度是否超過了最大飛翔速度
            if v(i,j)>vmax;
                v(i,j)=vmax;
            elseif v(i,j)<-vmax;
                v(i,j)=-vmax;
            end
        end
        x(i,:)=x(i,:)+v(i,:);
    end
    if abs(globalbest_faval)<E0,break,end
    k=k+1;
end
Value1=1/globalbest_faval-1; 
Value1=num2str(Value1);
% strcat指令可以實現字元的組合輸出
disp(strcat('the maximum value','=',Value1));
%輸出最大值所在的橫座標位置
Value2=globalbest_x; 
% % Value2=num2str(Value2);
% % disp(strcat('the corresponding coordinate','=',Value2));
% x=-2:0.1:3;
% y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);
% plot(x,y,'m-','linewidth',3);
% hold on;
% plot(globalbest_x,1/globalbest_faval-1,'kp','linewidth',4);
% legend('樣本的輸出值','網路的輸出值');xlabel('X');ylabel('Y');grid on;toc;

% 賦值給 bp神經網路的權值和閾值
W1=zeros(HiddenNum,InDim);
% B1=zeros(HiddenNum,1); 
% W2=zeros(1,HiddenNum);           
% B2=zeros(1,OutDim); 

    biaozhi=1;
    for wi=1:HiddenNum
        W1(wi,:)=Value2(1,biaozhi:biaozhi+InDim-1);
        biaozhi=wi*InDim+1;
    end
    B1=Value2(1,biaozhi:biaozhi+HiddenNum-1);
    B1=B1';
    W2=Value2(1,biaozhi+HiddenNum:biaozhi+HiddenNum+HiddenNum-1);
    B2=Value2(1,biaozhi+HiddenNum+HiddenNum);
end

 

相關文章