粒子群算法matlab源程序function

　　

%

function [pso F] = PSO_2D()

% FUNCTION PSO ---USE Particle Swarm Optimization Algorithm

%global present;

% close all;

pop_size = 10; % pop_size 種群大小

part_size = 2; % part_size 粒子大小, ** =n-D

gbest = zeros(1,part_size+1); % gbest 當前搜索到的最小的值

max_gen = 80; % max_gen 最大迭代次數

region=zeros(part_size,2); % 設定搜索空間範圍

region=[-3,3;-3,3]; % **每一維設定不同範圍

rand('state',sum(100*clock)); % 重置隨機數發生器狀態

arr_present = ini_pos(pop_size,part_size); % present 當前位置,隨機初始化,rand()的範圍為0~1

v=ini_v(pop_size,part_size); % 初始化當前速度

pbest = zeros(pop_size,part_size+1); % pbest 粒子以前搜索到的最優值，最後一列包括這些值的適應度

w_max = 0.9; % w_max 權係數最大值

w_min = 0.4;

v_max = 2; % **最大速度,為粒子的範圍寬度

c1 = 2; % 學習因子

c2 = 2; % 學習因子

best_record = zeros(1,max_gen); % best_record記錄最好的粒子的適應度。

% ————————————————————————

% 計算原始種群的適應度,及初始化

% ————————————————————————

arr_present(:,end)=ini_fit(arr_present,pop_size,part_size);

% for k=1:pop_size

% present(k,end) = fitness(present(k,1:part_size)); %計算原始種群的適應度

% end

pbest = arr_present; %初始化各個粒子最優值

[best_value best_index] = min(arr_present(:,end)); %初始化全局最優，即適應度為全局最小的值，根據需要也可以選取為最大值

gbest = arr_present(best_index,:);

%v = zeros(pop_size,1); % v 速度

% ————————————————————————

% 迭代

% ————————————————————————

% global m;

% m = moviein(1000); %生成幀矩陣

x=[-3:0.01:3];

y=[-3:0.01:3];

z=@(x,y) 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...

- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...

- 1/3*exp(-(x+1).^2 - y.^2);

for i=1:max_gen

grid on;

plot3(x,y,z);

% subplot(121),ezmesh(z),hold on,grid on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;

% subplot(122),ezmesh(z),view([145,90]),hold on,grid on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;

ezmesh(z),hold on,grid on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;

drawnow

F(i)=getframe;

% ezmesh(z)

% % view([-37,90])

% hold on;

% grid on;

% % plot(-0.0898,0.7126,'ro');

% plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'); %改為三維

% axis([-2*pi,2*pi,-pi,pi,-50,10]);

% hold off;

pause(0.01);

% m(:,i) = getframe; %添加圖形

w = w_max-(w_max-w_min)*i/max_gen;

% fprintf('# %i 代開始！\n',i);

% 確定是否對打散已經收斂的粒子群——————————————————————————————

reset = 0; % reset = 1時設置為粒子群過分收斂時將其打散，如果＝1則不打散

if reset==1

bit = 1;

for k=1:part_size

bit = bit&(range(arr_present(:,k))<0.1);

end

if bit==1 % bit=1時對粒子位置及速度進行隨機重置

arr_present = ini_pos(pop_size,part_size); % present 當前位置,隨機初始化

v = ini_v(pop_size,part_size); % 速度初始化

for k=1:pop_size % 重新計算適應度

arr_present(k,end) = fitness(arr_present(k,1:part_size));

end

warning('粒子過分集中！重新初始化……'); % 給出信息

display(i);

end

end

for j=1:pop_size

v(j,:) = w.*v(j,:)+c1.*rand.*(pbest(j,1:part_size)-arr_present(j,1:part_size))...

+c2.*rand.*(gbest(1:part_size)-arr_present(j,1:part_size)); % 粒子速度更新 (a)

% 判斷v的大小，限制v的絕對值小於5————————————————————————————

c = find(abs(v)>6); %**最大速度設置，粒子的範圍寬度

v(c) = sign(v(c))*6; %如果速度大於3.14則，速度為3.14

arr_present(j,1:part_size) = arr_present(j,1:part_size)+v(j,1:part_size); % 粒子位置更新 (b)

arr_present(j,end) = fitness(arr_present(j,1:part_size));

if (arr_present(j,end)>pbest(j,end))&(Region_in(arr_present(j,:),region)) % 根據條件更新pbest,如果是最小的值為小於號,相反則為大於號

pbest(j,:) = arr_present(j,:);

end

end

[best best_index] = max(arr_present(:,end)); % 如果是最小的值為min,相反則為max

if best>gbest(end)&(Region_in(arr_present(best_index,:),region)) % 如果當前最好的結果比以前的好，則更新最優值gbest,如果是最小的值為小於號,相反則為大於號

gbest = arr_present(best_index,:);

end

best_record(i) = gbest(end);

end

pso = gbest;

display(gbest);

% figure;

% plot(best_record);

% movie2avi(F,'pso_2D1.avi','compression','MSVC');

% ***************************************************************************

% 計算適應度

% ***************************************************************************

function fit = fitness(present)

fit=3*(1-present(1)).^2.*exp(-(present(1).^2) - (present(2)+1).^2) ... %**需要求極值的函數,本例即peaks函數

- 10*(present(1)/5 - present(1).^3 - present(2).^5).*exp(-present(1).^2-present(2).^2) ...

- 1/3*exp(-(present(1)+1).^2 - present(2).^2);

function ini_present=ini_pos(pop_size,part_size)

ini_present = 3*rand(pop_size,part_size+1); %初始化當前粒子位置,使其隨機的分布在工作空間 %** 6即為自變量範圍

function ini_velocity=ini_v(pop_size,part_size)

ini_velocity =3/2*(rand(pop_size,part_size)); %初始化當前粒子速度,使其隨機的分布在速度範圍內

function flag=Region_in(pos_present,region)

[m n]=size(pos_present);

flag=1;

for j=1:n-1

flag=flag&(pos_present(1:j)>=region(j,1))&(pos_present(1:j)<=region(j,2));

end

function arr_fitness=ini_fit(pos_present,pop_size,part_size)

for k=1:pop_size

arr_fitness(k,1) = fitness(pos_present(k,1:part_size)); %計算原始種群的適應度

end