spso.m

保证全局收敛的随机微粒群算法。当最优粒子的解无进化

function main()
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
POPSIZE = 20;
W = 0.73;
C1 = 2.05;
C2 = 2.05;
VMAX = 5;
% clusters spherical coordinates bounds
LJmin = [0,-1.0000,-3.0000,-6.0000,-9.1038,-12.7120,-16.5053,-19.8214,-24.1133,-28.4225,-32.7659,-37.9676,-44.3268];
% N = 3; radius = (3*N/(4*pi*sqrt(2)))^(1/3);
% unit = [0 2*pi;0 pi;0 radius];bounds=unit;
% for i=1:N-1 bounds = cat(1,bounds,unit); end
% clear radius unit
bounds = ones(2,1) * [-5.12 5.12];
XMIN = bounds(:,1)';
XMAX = bounds(:,2)';
DIMENSION = length(XMIN);
% DIMENSION = 2;
% XMIN = -5.12;
% XMAX = 5.12;
P=[]; % 1:DIMENSION
PBEST=[];
individual = struct([]);%struct('number','best','bestfitness','fitness','speed');

% fprintf(1,'Minimization of LJ clusters: N = %d\n',N);
total = 0;
sum = 0;
tic
for j=1:1
    initiate;
    for i=1:2000
        opt = PBEST;
        flag = -1;
        if(round(rand))
            direc = 1;
        else 
            direc = 1;
        end
        for k=1:POPSIZE
            if((individual(k).fitness == individual(k).bestfitness) && (individual(k).fitness == PBEST))
                flag = k;
                break;
            end
        end
        if(flag ~= -1)
            generate(flag);
        end       
        for k=1:POPSIZE
            if(k ~= flag)
                for t=1:DIMENSION
                    individual(k).speed(t) = W * ( individual(k).speed(t) + direc * ( C1 * rand * (individual(k).best(t) - individual(k).number(t)) + ...
                                                                                      C2 * rand * (P(t) - individual(k).number(t))) );
                    if (individual(k).speed(t) > VMAX)	individual(k).speed(t) = VMAX; end
                    individual(k).number(t) = individual(k).number(t) + individual(k).speed(t);
%                     if(individual(k).number(t) < XMIN(t)) 
%                         individual(k).number(t) = XMIN(t);
%                     elseif(individual(k).number(t) > XMAX(t)) 
%                         individual(k).number(t) = XMAX(t);
%                     end
                end
                calculation(k);
                localbest(k);
            end
%                     direc = calc_div(direc,i);
        end
        globalbest(1);
        if(abs(PBEST)<1e-3)
%         if(abs(PBEST-LJmin(N))<1e-3) 
            total = total + 1;
            sum = sum + i;
            break
        elseif((opt-PBEST)>1e-3)
            fprintf(1,'\n%d %.4f\n',i,PBEST);
        else
            fprintf(1,'%d ',i);           
        end
    end 
    fprintf(1,'\n%d\t%.4f\n',i,PBEST);    
end
fprintf(1,'Total number is %d\nSum is %d\nAverage is %.1f\n',total,sum,sum/total);
toc

function calculation(num)
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
x = individual(num).number;
% x1 = individual(num).number(2);%x0 = -1.3071;x1 = -1.4249;
% s = 0;
% h = 0;
% for i=1:5
%     s = s + i * cos((i-1)*x0+i);
%     h = h + i * cos((i+1)*x1+i);
% end
% individual(num).fitness = s * h + (x0+1.42513)^2.0 + (x1+0.80032)^2.0;
% s = (x0+x1+1.0)^2.0 * (19.0-14.0*x0+3.0*x0^2.0-14.0*x1+6.0*x0*x1+3.0*x1^2.0) + 1.0;
% h = 30.0 + (2.0*x0-3.0*x1)^2.0 * (18.0-32.0*x0+12.0*x0^2.0+48.0*x1-36.0*x0*x1+27.0*x1^2.0);
% individual(num).fitness = s * h;
% Rastrigin's function (-5.12 5.12)
len=length(x);A=8;
fval=0;
for i=1:len
    fval=fval+x(i)^2-A*cos(2*pi*x(i));
end
fval=fval+len*A;
% spherical coordinates
% x = individual(num).number;
% N = DIMENSION / 3;
% fval=0;
% for i=1:N-1
%     for j=i+1:N
%         r=(x(3*j)*cos(x(3*j-2))*cos(x(3*j-1))-x(3*i)*cos(x(3*i-2))*cos(x(3*i-1)))^2 + ...
%           (x(3*j)*cos(x(3*j-2))*sin(x(3*j-1))-x(3*i)*cos(x(3*i-2))*sin(x(3*i-1)))^2 + ...
%           (x(3*j)*sin(x(3*j-2))-x(3*i)*sin(x(3*i-2)))^2;
%         fval = fval + r^(-6) - 2*r^(-3);
%     end
% end
individual(num).fitness = fval;
return

function generate(num)
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
individual(num).number = rand(1,DIMENSION).*(XMAX-XMIN) + XMIN;
calculation(num);
individual(num).bestfitness = individual(num).fitness;
individual(num).best = individual(num).number;
return

function globalbest(number)
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
if(number == 0)
    [PBEST,flag] = max([individual.fitness]);
    P = individual(flag).number;
else
    for i=1:POPSIZE
        if(individual(i).bestfitness < PBEST)
            P = individual(i).best;
            PBEST = individual(i).bestfitness;
        end
    end
end
return

function initiate()
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
for i=1:POPSIZE
    individual(i).number = rand(1,DIMENSION).*(XMAX-XMIN) + XMIN;
    individual(i).speed = rand(1,DIMENSION).* VMAX;
    individual(i).best = individual(i).number;
    calculation(i);
    individual(i).bestfitness = individual(i).fitness;
end
globalbest(0);
return

function localbest(num)
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
if(individual(num).bestfitness > individual(num).fitness)
    individual(num).best = individual(num).number;
end
individual(num).bestfitness = individual(num).fitness;
return

function Stretch_Tech
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
gma1=1e+4/2;gma2=1/2;meu=1e-10;
for i=1:POPSIZE
    if(abs(individual(i).fitness-PBEST)>1e-5)
        x = individual(i).number;
        fval = individual(i).fitness;
        G = fval + gma1 * sqrt(sum((x-P).^2)) * (sign(fval-PBEST)+1);
        H = G + gma2 * (sign(fval-PBEST)+1) / tanh(meu*(G-PBEST));
        individual(i).fitness = H;
    end
end
return

function direc = calc_div(direc,iter)
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
D_low = 5e-6 * exp(1-iter/2000);
D_high = 0.25 * exp(1-iter/2000);
div = 0;
x = [individual.number];
x = reshape(x,DIMENSION,POPSIZE);
x = x';
for i=1:POPSIZE
    for j=1:DIMENSION
        div = div + (individual(i).number(j) - mean(x(:,j)))^2;
    end
end
div = div / ( POPSIZE * max(XMAX-XMIN) );
if(direc>0 && div<D_low)
    direc = -1;
elseif(direc<0 && div>D_high)
    direc = 1;
end
return

function direc = calc_dir2(t)
global POPSIZE DIMENSION W C1 C2 VMAX XMIN XMAX P PBEST individual
D_low = 5e-6;
D_high = 0.25;
div = 0;
x = [individual.number];
x = reshape(x,DIMENSION,POPSIZE);
x = x';
for i=1:POPSIZE
    div = div + (individual(i).number(t) - mean(x(:,t)))^2;
end
div = div / (POPSIZE * (XMAX(t)-XMIN(t)));
if(div<D_low)
    direc = -1;
elseif(div>D_high)
    direc = 1;
elseif round(rand) 
    direc = -1;
else 
    direc = 1;
end
return