​

一、粒子群演算法

粒子群演算法是在1995年由Eberhart博士和Kennedy博士一起提出的，它源於對鳥群捕食行為的研究。它的基本核心是利用群體中的個體對資訊的共享從而使整個群體的運動在問題求解空間中產生從無序到有序的演化過程，從而獲得問題的最優解。設想這麼一個場景：一群鳥進行覓食，而遠處有一片玉米地，所有的鳥都不知道玉米地到底在哪裡，但是它們知道自己當前的位置距離玉米地有多遠。那麼找到玉米地的最佳策略，也是最簡單有效的策略就是搜尋目前距離玉米地最近的鳥群的周圍區域。

在PSO中，每個優化問題的解都是搜尋空間中的一隻鳥，稱之為"粒子"，而問題的最優解就對應於鳥群中尋找的"玉米地"。所有的粒子都具有一個位置向量（粒子在解空間的位置）和速度向量（決定下次飛行的方向和速度），並可以根據目標函式來計算當前的所在位置的適應值（fitness value）,可以將其理解為距離"玉米地"的距離。在每次的迭代中，種群中的例子除了根據自身的經驗（歷史位置）進行學習以外，還可以根據種群中最優粒子的"經驗"來學習，從而確定下一次迭代時需要如何調整和改變飛行的方向和速度。就這樣逐步迭代，最終整個種群的例子就會逐步趨於最優解。

上面的解釋可能還比較抽象，下面通過一個簡單的例子來進行說明

​

在一個湖中有兩個人他們之間可以通訊，並且可以探測到自己所在位置的最低點。初始位置如上圖所示，由於右邊比較深，因此左邊的人會往右邊移動一下小船。

​

現在左邊比較深，因此右邊的人會往左邊移動一下小船

一直重複該過程，最後兩個小船會相遇

​

得到一個區域性的最優解

​

​

將每個個體表示為粒子。每個個體在某一時刻的位置表示為，x(t),方向表示為v(t)

​

p（t）為在t時刻x個體的自己的最優解，g(t)為在t時刻所有個體的最優解，v(t)為個體在t時刻的方向，x(t)為個體在t時刻的位置

​​

下一個位置為上圖所示由x,p,g共同決定了

​

種群中的粒子通過不斷地向自身和種群的歷史資訊進行學習，從而可以找到問題的最優解。

但是，在後續的研究中表表明，上述原始的公式中存在一個問題：公式中V的更新太具有隨機性，從而使整個PSO演算法的全域性優化能力很強，但是區域性搜尋能力較差。而實際上，我們需要在演算法迭代初期PSO有著較強的全域性優化能力，而在演算法的後期，整個種群應該具有更強的區域性搜尋能力。所以根據上述的弊端，shi和Eberhart通過引入慣性權重修改了公式，從而提出了PSO的慣性權重模型：

每一個向量的分量表示如下

​

​

其中w稱為是PSO的慣性權重，它的取值介於【0,1】區間，一般應用中均採用自適應的取值方法，即一開始令w=0.9，使得PSO全域性優化能力較強，隨著迭代的深入，引數w進行遞減，從而使的PSO具有較強的區域性優化能力，當迭代結束時，w=0.1。引數c1和c2稱為學習因子，一般設定為1,4961；而r1和r2為介於[0,1]之間的隨機概率值。

整個粒子群優化演算法的演算法框架如下：

step1種群初始化，可以進行隨機初始化或者根據被優化的問題設計特定的初始化方法，然後計算個體的適應值，從而選擇出個體的區域性最優位置向量和種群的全域性最優位置向量。

step2 迭代設定：設定迭代次數，並令當前迭代次數為1

step3 速度更新：更新每個個體的速度向量

step4 位置更新：更新每個個體的位置向量

step5 區域性位置和全域性位置向量更新：更新每個個體的區域性最優解和種群的全域性最優解

step6 終止條件判斷：判斷迭代次數時都達到最大迭代次數，如果滿足，輸出全域性最優解，否則繼續進行迭代，跳轉至step 3。

　　對於粒子群優化演算法的運用，主要是對速度和位置向量迭代運算元的設計。迭代運算元是否有效將決定整個PSO演算法效能的優劣，所以如何設計PSO的迭代運算元是PSO演算法應用的研究重點和難點。

二、競爭學習的粒子群優化演算法(CLPSO)

傳統的時變慣性權重粒子群優化演算法對於求解一般的全域性最優問題具有良好的效果,而對於複雜高維的優化問題易陷入區域性收斂、存在早熟等缺點。針對以上存在的缺點,提出了種群進化離散度的概念,並考慮Sigmoid函式線上性與非線性之間較好的平衡效能,給出一種非線性動態自適應慣性權重的粒子群優化演算法。該演算法充分考慮進化過程中種群粒子之間進化差異,自適應地賦予不同的慣性權重因子,滿足粒子群優化演算法在不同進化時期對全域性探索和區域性開發能力的需求,模擬例項測試結果驗證了該演算法的有效性。

三、部分程式碼

function [gbest,gbestval,fitcount]= CLPSO_new_func(fhd,Max_Gen,Max_FES,Particle_Number,Dimension,VRmin,VRmax,varargin)
%[gbest,gbestval,fitcount]= CLPSO_new_func('f8',3500,200000,30,30,-5.12,5.12)
rand('state',sum(100*clock));
me=Max_Gen;
ps=Particle_Number;
D=Dimension;
cc=[1 1];   %acceleration constants
t=0:1/(ps-1):1;t=5.*t;
Pc=0.0+(0.5-0.0).*(exp(t)-exp(t(1)))./(exp(t(ps))-exp(t(1)));
% Pc=0.5.*ones(1,ps);
m=0.*ones(ps,1);
iwt=0.9-(1:me)*(0.7/me);
% iwt=0.729-(1:me)*(0.0/me);
cc=[1.49445 1.49445];
if length(VRmin)==1
    VRmin=repmat(VRmin,1,D);
    VRmax=repmat(VRmax,1,D);
end
mv=0.2*(VRmax-VRmin);
VRmin=repmat(VRmin,ps,1);
VRmax=repmat(VRmax,ps,1);
Vmin=repmat(-mv,ps,1);
Vmax=-Vmin;
pos=VRmin+(VRmax-VRmin).*rand(ps,D);

for i=1:ps;
    e(i,1)=feval(fhd,pos(i,:),varargin{:});
end

fitcount=ps;
vel=Vmin+2.*Vmax.*rand(ps,D);%initialize the velocity of the particles
pbest=pos;
pbestval=e; %initialize the pbest and the pbest's fitness value
[gbestval,gbestid]=min(pbestval);
gbest=pbest(gbestid,:);%initialize the gbest and the gbest's fitness value
gbestrep=repmat(gbest,ps,1);

stay_num=zeros(ps,1);

ai=zeros(ps,D);
f_pbest=1:ps;f_pbest=repmat(f_pbest',1,D);
for k=1:ps
    ar=randperm(D);
    ai(k,ar(1:m(k)))=1;
    fi1=ceil(ps*rand(1,D));
    fi2=ceil(ps*rand(1,D));
    fi=(pbestval(fi1)<pbestval(fi2))'.*fi1+(pbestval(fi1)>=pbestval(fi2))'.*fi2;
    bi=ceil(rand(1,D)-1+Pc(k));
    if bi==zeros(1,D),rc=randperm(D);bi(rc(1))=1;end
    f_pbest(k,:)=bi.*fi+(1-bi).*f_pbest(k,:);
end

stop_num=0;
i=1;


while i<=me&fitcount<=Max_FES
    i=i+1;
    for k=1:ps
        
        if stay_num(k)>=5
            %     if round(i/10)==i/10%|stay_num(k)>=5
            stay_num(k)=0;
            ai(k,:)=zeros(1,D);
            f_pbest(k,:)=k.*ones(1,D);
            ar=randperm(D);
            ai(k,ar(1:m(k)))=1;
            fi1=ceil(ps*rand(1,D));
            fi2=ceil(ps*rand(1,D));
            fi=(pbestval(fi1)<pbestval(fi2))'.*fi1+(pbestval(fi1)>=pbestval(fi2))'.*fi2;
            bi=ceil(rand(1,D)-1+Pc(k));
            if bi==zeros(1,D),rc=randperm(D);bi(rc(1))=1;end
            f_pbest(k,:)=bi.*fi+(1-bi).*f_pbest(k,:);
        end
        
        for dimcnt=1:D
            pbest_f(k,dimcnt)=pbest(f_pbest(k,dimcnt),dimcnt);
        end
        aa(k,:)=cc(1).*(1-ai(k,:)).*rand(1,D).*(pbest_f(k,:)-pos(k,:))+cc(2).*ai(k,:).*rand(1,D).*(gbestrep(k,:)-pos(k,:));%~~~~~~~~~~~~~~~~~~~~~~
        vel(k,:)=iwt(i).*vel(k,:)+aa(k,:);
        vel(k,:)=(vel(k,:)>mv).*mv+(vel(k,:)<=mv).*vel(k,:);
        vel(k,:)=(vel(k,:)<(-mv)).*(-mv)+(vel(k,:)>=(-mv)).*vel(k,:);
        pos(k,:)=pos(k,:)+vel(k,:);
        
        if (sum(pos(k,:)>VRmax(k,:))+sum(pos(k,:)<VRmin(k,:)))==0;
            e(k,1)=feval(fhd,pos(k,:),varargin{:});
            fitcount=fitcount+1;
            tmp=(pbestval(k)<=e(k));
            if tmp==1
                stay_num(k)=stay_num(k)+1;
            end
            temp=repmat(tmp,1,D);
            pbest(k,:)=temp.*pbest(k,:)+(1-temp).*pos(k,:);
            pbestval(k)=tmp.*pbestval(k)+(1-tmp).*e(k);%update the pbest
            if pbestval(k)<gbestval
                gbest=pbest(k,:);
                gbestval=pbestval(k);
                gbestrep=repmat(gbest,ps,1);%update the gbest
            end
        end
        
    end
    
    % if round(i/100)==i/100
    %     plot(pos(:,D-1),pos(:,D),'b*');hold on;
    %     for k=1:floor(D/2)
    %         plot(gbest(:,2*k-1),gbest(:,2*k),'r*');
    %     end
    %     hold off
    %     title(['PSO: ',num2str(i),' generations, Gbestval=',num2str(gbestval)]);
    %     axis([VRmin(1,D-1),VRmax(1,D-1),VRmin(1,D),VRmax(1,D)])
    %     drawnow
    % end
    
    if fitcount>=Max_FES
        break;
    end
    if (i==me)&(fitcount<Max_FES)
        i=i-1;
    end
end
gbestval

四、模擬結果

​

五、參考文獻及程式碼私信博主

[1]王生亮,劉根友.一種非線性動態自適應慣性權重PSO演算法[J].計算機模擬,2021,38(04):249-253+451.

​