【優化演算法】粒子群工具箱函式優化演算法【含Matlab原始碼 1126期】
一、簡介
粒子群演算法源於複雜適應系統(Complex Adaptive System,CAS)。CAS理論於1994年正式提出,CAS中的成員稱為主體。比如研究鳥群系統,每個鳥在這個系統中就稱為主體。主體有適應性,它能夠與環境及其他的主體進行交流,並且根據交流的過程“學習”或“積累經驗”改變自身結構與行為。整個系統的演變或進化包括:新層次的產生(小鳥的出生);分化和多樣性的出現(鳥群中的鳥分成許多小的群);新的主題的出現(鳥尋找食物過程中,不斷髮現新的食物)。
所以CAS系統中的主體具有4個基本特點(這些特點是粒子群演算法發展變化的依據):
首先,主體是主動的、活動的。
主體與環境及其他主體是相互影響、相互作用的,這種影響是系統發展變化的主要動力。
環境的影響是巨集觀的,主體之間的影響是微觀的,巨集觀與微觀要有機結合。
最後,整個系統可能還要受一些隨機因素的影響。
粒子群演算法就是對一個CAS系統---鳥群社會系統的研究得出的。
粒子群演算法( Particle Swarm Optimization, PSO)最早是由Eberhart和Kennedy於1995年提出,它的基本概念源於對鳥群覓食行為的研究。設想這樣一個場景:一群鳥在隨機搜尋食物,在這個區域裡只有一塊食物,所有的鳥都不知道食物在哪裡,但是它們知道當前的位置離食物還有多遠。那麼找到食物的最優策略是什麼呢?最簡單有效的就是搜尋目前離食物最近的鳥的周圍區域。
PSO演算法就從這種生物種群行為特性中得到啟發並用於求解優化問題。在PSO中,每個優化問題的潛在解都可以想象成d維搜尋空間上的一個點,我們稱之為“粒子”(Particle),所有的粒子都有一個被目標函式決定的適應值(Fitness Value ),每個粒子還有一個速度決定他們飛翔的方向和距離,然後粒子們就追隨當前的最優粒子在解空間中搜索。Reynolds對鳥群飛行的研究發現。鳥僅僅是追蹤它有限數量的鄰居但最終的整體結果是整個鳥群好像在一箇中心的控制之下.即複雜的全域性行為是由簡單規則的相互作用引起的。
二、原始碼
%% 基於粒子群工具箱的函式優化演算法 %% 清空環境 clear clc %% 引數初始化 x_range=[-50,50]; %引數x變化範圍 y_range=[-50,50]; %引數y變化範圍 range = [x_range;y_range]; %引數變化範圍(組成矩陣) Max_V = 0.2*(range(:,2)-range(:,1)); %最大速度取變化範圍的10%~20% n=2; %待優化函式的維數,此例子中僅x、y兩個自變數,故為2 PSOparams= [25 2000 24 2 2 0.9 0.4 1500 1e-25 250 NaN 0 0]; %% 粒子群尋優 pso_Trelea_vectorized('Rosenbrock',n,Max_V,range,0,PSOparams) %呼叫PSO核心模組 % goplotpso.m % default plotting script used in PSO functions % % this script is not a function, % it is a plugin for the main PSO routine (pso_Trelea_vectorized) % so it shares all the same variables, be careful with variable names % when making your own plugin % setup figure, change this for your own machine clf set(gcf,'Position',[651 31 626 474]); % this is the computer dependent part %set(gcf,'Position',[743 33 853 492]); set(gcf,'Doublebuffer','on'); % particle plot, upper right subplot('position',[.7,.6,.27,.32]); set(gcf,'color','k') plot3(pos(:,1),pos(:,D),out,'b.','Markersize',7) hold on plot3(pbest(:,1),pbest(:,D),pbestval,'g.','Markersize',7); plot3(gbest(1),gbest(D),gbestval,'r.','Markersize',25); % crosshairs offx = max(abs(min(min(pbest(:,1)),min(pos(:,1)))),... abs(max(max(pbest(:,1)),max(pos(:,1))))); offy = max(abs(min(min(pbest(:,D)),min(pos(:,D)))),... abs(min(max(pbest(:,D)),max(pos(:,D))))); plot3([gbest(1)-offx;gbest(1)+offx],... [gbest(D);gbest(D)],... [gbestval;gbestval],... 'r-.'); plot3([gbest(1);gbest(1)],... [gbest(D)-offy;gbest(D)+offy],... [gbestval;gbestval],... 'r-.'); hold off xlabel('Dimension 1','color','y') ylabel(['Dimension ',num2str(D)],'color','y') zlabel('Cost','color','y') title('Particle Dynamics','color','w','fontweight','bold') set(gca,'Xcolor','y') set(gca,'Ycolor','y') set(gca,'Zcolor','y') set(gca,'color','k') % camera control view(2) try axis([gbest(1)-offx,gbest(1)+offx,gbest(D)-offy,gbest(D)+offy]); catch axis([VR(1,1),VR(1,2),VR(D,1),VR(D,2)]); end % error plot, left side subplot('position',[0.1,0.1,.475,.825]); semilogy(tr(find(~isnan(tr))),'color','m','linewidth',2) %plot(tr(find(~isnan(tr))),'color','m','linewidth',2) xlabel('epoch','color','y') ylabel('gbest val.','color','y') if D==1 titstr1=sprintf(['%11.6g = %s( [ %9.6g ] )'],... gbestval,strrep(functname,'_','\_'),gbest(1)); elseif D==2 titstr1=sprintf(['%11.6g = %s( [ %9.6g, %9.6g ] )'],... gbestval,strrep(functname,'_','\_'),gbest(1),gbest(2)); elseif D==3 titstr1=sprintf(['%11.6g = %s( [ %9.6g, %9.6g, %9.6g ] )'],... gbestval,strrep(functname,'_','\_'),gbest(1),gbest(2),gbest(3)); else titstr1=sprintf(['%11.6g = %s( [ %g inputs ] )'],... gbestval,strrep(functname,'_','\_'),D); end title(titstr1,'color','m','fontweight','bold'); grid on % axis tight set(gca,'Xcolor','y') set(gca,'Ycolor','y') set(gca,'Zcolor','y') set(gca,'color','k') set(gca,'YMinorGrid','off') % text box in lower right % doing it this way so I can format each line any way I want subplot('position',[.62,.1,.29,.4]); clear titstr if trelea==0 PSOtype = 'Common PSO'; xtraname = 'Inertia Weight : '; xtraval = num2str(iwt(length(iwt))); elseif trelea==2 | trelea==1 PSOtype = (['Trelea Type ',num2str(trelea)]); xtraname = ' '; xtraval = ' '; elseif trelea==3 PSOtype = (['Clerc Type 1"']); xtraname = '\chi value : '; xtraval = num2str(chi); end if isnan(errgoal) errgoalstr='Unconstrained'; else errgoalstr=num2str(errgoal); end if minmax==1 minmaxstr = ['Maximize to : ']; elseif minmax==0 minmaxstr = ['Minimize to : ']; else minmaxstr = ['Target to : ']; end if rstflg==1 rststat1 = 'Environment Change'; rststat2 = ' '; else rststat1 = ' '; rststat2 = ' '; end titstr={'PSO Model: ' ,PSOtype;... 'Dimensions : ' ,num2str(D);... '# of particles : ',num2str(ps);... minmaxstr ,errgoalstr;... 'Function : ' ,strrep(functname,'_','\_');... xtraname ,xtraval;... rststat1 ,rststat2}; text(.1,1,[titstr{1,1},titstr{1,2}],'color','g','fontweight','bold'); hold on text(.1,.9,[titstr{2,1},titstr{2,2}],'color','m'); text(.1,.8,[titstr{3,1},titstr{3,2}],'color','m'); text(.1,.7,[titstr{4,1}],'color','w'); text(.55,.7,[titstr{4,2}],'color','m'); text(.1,.6,[titstr{5,1},titstr{5,2}],'color','m'); text(.1,.5,[titstr{6,1},titstr{6,2}],'color','w','fontweight','bold'); text(.1,.4,[titstr{7,1},titstr{7,2}],'color','r','fontweight','bold'); % if we are training a neural net, show a few more parameters if strcmp('pso_neteval',functname) % net is passed from trainpso to pso_Trelea_vectorized in case you are % wondering where that structure comes from hiddlyrstr = []; for lyrcnt=1:length(net.layers) TF{lyrcnt} = net.layers{lyrcnt}.transferFcn; Sn(lyrcnt) = net.layers{lyrcnt}.dimensions; hiddlyrstr = [hiddlyrstr,', ',TF{lyrcnt}]; end hiddlyrstr = hiddlyrstr(3:end); text(0.1,.35,['#neur/lyr = [ ',num2str(net.inputs{1}.size),' ',... num2str(Sn),' ]'],'color','c','fontweight','normal',... 'fontsize',10); text(0.1,.275,['Lyr Fcn: ',hiddlyrstr],... 'color','c','fontweight','normal','fontsize',9); end legstr = {'Green = Personal Bests';... 'Blue = Current Positions';... 'Red = Global Best'}; text(.1,0.025,legstr{1},'color','g'); text(.1,-.05,legstr{2},'color','b'); text(.1,-.125,legstr{3},'color','r'); hold off set(gca,'color','k'); set(gca,'visible','off'); drawnow
三、執行結果
四、備註
版本:2014a
完整程式碼或代寫加QQ1564658423