基于PSO三維極點搜尋源碼

function [pso F] = pso()

%FUNCTION Optimization --------USE Particle Swarm Optimization Algorithm

% % %

% Msc. Thesis

% Rentian Huang

% January, 2007

%global present;

% close all;

pop_size = 10;          %   pop_size 

part_size = 2;          %   partible_size ,   ** =n-D

gbest = zeros(1,part_size+1);            %   gbest the best value so far

max_gen = 80;          %   max_gen max generation

region=zeros(part_size,2);  %   define searching regions

region=[-1,1;-1,1];          %    define the size of regions

rand('state',sum(100*clock));   %   reset the state of random machine

arr_present = ini_pos(pop_size,part_size);   %   present position, random initialization,rand() from 0~1

v=ini_v(pop_size,part_size);             %   initialized velocity

pbest = zeros(pop_size,part_size+1);      %   pbest the best result, last column contain the parameter as well

w_max = 0.9;                            %   w_max maximum of weight

w_min = 0.4;                            %   w_min minumum of weitht

v_max = 2;             %      ** maximum of velocity 

c1 = 2;                   %   learning parameter 1

c2 = 2;                   %   learning parameter 2

best_record = zeros(1,max_gen);     %   best_record remember the best particle's parameters

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

%   Calulate the paricles affinity and initialization

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

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;                                        %initialized all best value for particles

[best_value best_index] = min(arr_present(:,end));         %initialize the global optima

gbest = arr_present(best_index,:);

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

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

%   Evolving

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

% global m;

% m = moviein(1000);      %

x=[-1:0.01:1];

y=[-1:0.01:1];

[email protected](x,y) cos(2 * pi .* x) .* cos(2 * pi .* y) .* exp(-0.1 * (x.^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,[-1,1,-1,1]),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);

    %   Decide whether reset the particles because they may converge——————————————————————————————

    reset = 0;          %   reset = 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, reset position and velocity

            arr_present = ini_pos(pop_size,part_size);   %   present position

            v = ini_v(pop_size,part_size);           %   velocity initialization

            for k=1:pop_size                                    %   re-calculate affinity

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

            end

            warning('Reset particle because they already converged……');    

            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));                        %  renew velocity (a)

        %   The abs(velocity) must <5————————————————————————————

        c = find(abs(v)>6);                                                                 

        v(c) = sign(v(c))*6;   %if v>3.14則，v=3.14

        arr_present(j,1:part_size) = arr_present(j,1:part_size)+v(j,1:part_size);              %  renew position (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))     %   renew pbest

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

        end

    end

    [best best_index] = max(arr_present(:,end));                                     

    if best>gbest(end)&(Region_in(arr_present(best_index,:),region))   %   if the result better than before, save to 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');