apso.m

改进的粒子群算法--自适应粒子群算法

function [pso t i] = apso(max_gen,pop_size,part_size)
%   FUNCTION PSO  --------USE Particle Swarm Optimization Algorithm

%clear all
t=cputime;
D_max = 0.25;
D_min = 0.001;
E_max = 2;
E_min = 0.25;
Pm = 0.005;
%pop_size = 30;          %   pop_size 种群大小
%part_size = 30;          %   part_size 粒子大小, 
Q  =  10*part_size;                   %   划分的块数
gbest = zeros(1,part_size+1);            %   gbest 当前搜索到的最小的值
%max_gen = 30;          %   max_gen 最大迭代次数
region=zeros(part_size,2);  %   设定搜索空间范围
%region=[-300,300;-300,300;-300,300;-300,300;-300,300;-300,300;-300,300;-300,300;-300,300;-300,300];          %区域                                                                                      **每一维设定不同范围
% region=[-30,30;-30,30;-30,30;-30,30;-30,30;-30,30;-30,30;-30,30;-30,30;-30,30];  
 %region=[-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3];
 region=[-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3];
mode='exploit';
rand('state',sum(100*clock));   %   重置随机数发生器状态
arr_present = ini_pos(pop_size,part_size,region);   %   present 当前位置,随机初始化,rand()的范围为0~1
v_max = 1;   
v=ini_v(pop_size,part_size,v_max);             %   初始化当前速度


pbest = zeros(pop_size,part_size+1);      %   pbest 粒子以前搜索到的最优值，最后一列包括这些值的适应度
w_max = 0.7;                            %   w_max 权系数最大值
w_min = 0.05;
          %                                                                                          **最大速度,为粒子的范围宽度
c1 = 2;                   %   学习因子
c2 = 2;                   %   学习因子
best_record = zeros(1,max_gen);     %   best_record记录最好的粒子的适应度。
%  ————————————————————————
%   计算原始种群的适应度,及初始化
%  ————————————————————————
arr_present(:,end)=ini_fit(arr_present,pop_size,part_size);
pbest = arr_present;                                        %初始化各个粒子最优值
[best_value best_index] = min(arr_present(:,end));         %初始化全局最优，即适应度为全局最小的值，根据需要也可以选取为最大值
gbest = arr_present(best_index,:);
% i=1;
% while (i<=max_gen&abs(gbest(end))>0.01)
%   i=i+1;
  for i=1:max_gen  
    D = D_caculation(arr_present,region);
    E = E_caculation(arr_present,Q,pop_size,region);
    if D < D_min || E<E_min
        mode='explore';
    end
    if D > D_max && E>E_max
        mode='exploit';
    end
    

   
   if mode == 'exploit'
       w = D*(w_max-w_min)/(D_max-D_min)+(D_max*w_min-D_min*w_max)/(D_max-D_min);
      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的绝对值小于v_max————————————————————————————
        c = find(abs(v)>v_max);                                                                                              %**最大速度设置，粒子的范围宽度
        v(c) = sign(v(c))*v_max;   

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

        if (arr_present(j,end)<pbest(j,end))     %   根据条件更新pbest,如果是最小的值为小于号,相反则为大于号
           pbest(j,:) = arr_present(j,:);
           
        end

      end
   else 
       for j=1:pop_size
          if rand<Pm 
             arr_present(j,1:part_size+1)=(region(1,2)-region(1,1))*(rand(1,part_size+1)-0.5*ones(1,part_size+1)); 
             v(j,:)=v_max*(rand(1,part_size)-0.5*ones(1,part_size));
         
              arr_present(j,end) = fitness(arr_present(j,1:part_size));
             if (arr_present(j,end)<pbest(j,end))   %   根据条件更新pbest,如果是最小的值为小于号,相反则为大于号
               pbest(j,:) = arr_present(j,:);
             end
           
          end
          
       end
      
   end
    

    [best best_index] = min(pbest(:,end));                                      %   如果是最小的值为min,相反则为max

    if best<gbest(end)                %   如果当前最好的结果比以前的好，则更新最优值gbest,如果是最小的值为小于号,相反则为大于号
        gbest = pbest(best_index,:);
    end

    best_record(i) = gbest(end);
    
end
% display(i)
pso = gbest;

display(gbest);

t=cputime-t



%   ***************************************************************************
%      计算适应度
%   ***************************************************************************
function fit = fitness(present)
[m n]=size(present);
fit=0;
temp=1;
for ik=1:n
    fit=fit+present(:,ik)^2;
    temp=temp*cos(present(:,ik)/ik);
end
fit=fit/4000-temp+1;

% [m n]=size(present);
% fit=20*exp(-0.2*sqrt(sum(present.^2)/n))+exp(1/n*sum(cos(present.*2*3.1415926)))-20-2.72828;
%%%%%%%%%%%
% 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,region)
ini_present = (region(1,2)-region(1,1))*(rand(pop_size,part_size+1)-0.5*ones(pop_size,part_size+1));        %初始化当前粒子位置,使其随机的分布在工作空间                         %** 6即为自变量范围

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ini_velocity=ini_v(pop_size,part_size,v_max)
ini_velocity =v_max*(rand(pop_size,part_size)-0.5*ones(pop_size,part_size));        %初始化当前粒子速度,使其随机的分布在速度范围内


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pos_present=constr_pos(pos_present,region,part_size)
         for ik = 1 : part_size
             if pos_present(1,ik)<region(ik,1)
                 pos_present(1,ik)=region(ik,1);
             end
              if pos_present(1,ik)>region(ik,2)
                 pos_present(1,ik)=region(ik,2);
             end
             
         end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
function E = E_caculation(pos_position,Q,pop_size,region) 
Z=zeros(Q);

    for ik = 1:pop_size
         Z(floor((pos_position(ik,1)-region(1,1))/(region(1,2)-region(1,1))*Q)+1)=Z(floor((pos_position(ik,1)-region(1,1))/(region(1,2)-region(1,1))*Q)+1)+1;

     end
Z=Z/pop_size;
c=find(Z~=0);
E=-sum(Z(c).*log(Z(c)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = D_caculation(pos_position,region)
[pop_size part_size] = size(pos_position);
part_size = part_size-1;

for ikk = 1:part_size

      
    pos_avr_D(1,ikk) = sum(pos_position(:,ikk))/pop_size;
end

    for ik = 1:pop_size
       pos_D(ik,1) = sqrt(sum((pos_position(ik,1:part_size)-pos_avr_D(1,:)).^2));
    end
  D = sum(pos_D(:,1))/pop_size/(2*sqrt(sum(region(:,2).^2)));