【優化求解】基於慣性權值非線性遞減的改進粒子群演算法(IMPSO)matlab原始碼
一、粒子群演算法
粒子群演算法是在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演算法應用的研究重點和難點。
二、基於慣性權值非線性遞減的改進粒子群演算法(IMPSO)
針對粒子群優化演算法中出現的收斂早熟和不收斂的問題,提出了一種基於自然選擇和慣性權值非線性遞減的改進粒子群演算法,在演算法迭代過程中,粒子邊界速度採用最大速度非線性遞減變化策略來限制,慣性權值非線性遞減變化用於平衡種群粒子前期全域性搜尋與後期區域性尋優的能力;為使種群在進化過程中保持多樣性,在標準粒子群演算法中引用二階振盪策略使種群在進化過程中始終保持著多樣性;在此基礎上,進一步地將遺傳演算法中的選擇機理與粒子群演算法結合起來用於提高演算法的適用性能;所提出的演算法經過多個基準測試函式的模擬實驗驗證,並與其他已有演算法進行了對比;實驗結果表明:演算法在搜尋精度與尋優能力上有更明顯的優勢,尤其是在多維、多峰等複雜非線性優化問題時,所提演算法具有很強的競爭力。
三、部分程式碼
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] 華勇,王雙園,白國振,李炳初. 基於慣性權值非線性遞減的改進粒子群演算法 , 重慶工商大學學報(自然科學版), 2021 (02).