so_pso_xuan.m

一种带有灰色关联方法的粒子群算法

function [pso F] = pso_2D()
% FUNCTION PSO  --------USE Particle Swarm Optimization Algorithm
% global present;
% close all;
clc;
clear all;
pop_size = 30;                  %   pop_size 种群大小 ///粒子数量
part_size = 5;                  %   part_size 粒子大小 ///粒子的维数？
gbest = zeros(1,part_size+1);   %   gbest 当前搜索到的最小的值 
max_gen = 2000;                 %   max_gen 最大迭代次数

region=zeros(part_size,2);      % 设定搜索空间范围->解空间
region=100*[-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3]; % 每一维设定不同范围(称之为解空间，不是可行域空间）

rand('state',sum(100*clock));   % 重置随机数发生器状态


%当前种群的信息矩阵，逐代进化的群体 % 当前位置,随机初始化
% 一个30*10+1的随机的矩阵（初始化所有粒子的所有维数的位置值），其中最后一列为  
arr_present = ini_pos(pop_size,part_size); 
          
% 初始化当前速度 
% 一个10*2的随机的矩阵（初始化所有粒子的所有维数的速度值） 
v=ini_v(pop_size,part_size);  

%不是当前种群，可看作是一个外部的记忆体，存储每个粒子历史最优值（2维数值）：根据适应度更新！
%注意：pbest数组10*3  最后一列保存的是适应度
pbest = zeros(pop_size,part_size+1);    % pbest：粒子以前搜索到的最优值，最后一列包括这些值的适应度

% 1*80 保存每代的最优值
best_record = zeros(1,max_gen);         % best_record数组：记录每一代的最好的粒子的适应度


w_max = 0.9;          %  w_max权系数最大值
w_min = 0.4;          %  w_min权系数最小值
v_max = 2;            %  最大速度,为粒子的范围宽度
c1 = 2;               %  学习因子1
c2 = 2;               %  学习因子2

% ————————————————————————
%   计算原始种群的适应度,及初始化
% ————————————————————————
% 注意：传入的第一个参数是当前的粒子群体 ，ini_fit函数计算每个粒子的适应度
% arr_present(:,end)是最后一列 ，保存每个粒子的适应值，是这样的！xuan
arr_present(:,end)= ini_fit( arr_present, pop_size, part_size );

% 数组赋值，初始化每个粒子个体的历史最优值，以后会更新的 
pbest = arr_present;         % 初始化各个粒子最优值

% 找到当前群体中适应度最小的（在最后一列中寻找），best_value
% 改为max，表示关联度最大
[best_value best_index] = min(arr_present(:,end)); %初始化全局最优，即适应度为全局最小的值，根据需要也可以选取为最大值

% 唯一的全局最优值，是当前代所有粒子中最好的一个
gbest = arr_present(best_index,:);


% 因为是多目标，因此这个-----------------
%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);

%z1=(-10)*exp((-0.2)*sqrt(x(i)^2+x(i+1)^2));
%z2=(abs(x(i)))^0.8+5*sin(x(i)^3);


%开始进化，直到最大代数截至
for i=1:max_gen
    %grid on;
    %三维图象 %多维图象是画不出来的
    %ezmesh(z),hold on,grid on;
    %画出粒子群
    %plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;
    %drawnow 
    %flush
    %pause(0.01);
    
    w = w_max-(w_max-w_min)*i/max_gen; % 线形递减权重

    % 确定是否对打散已经收敛的粒子群—————————————————————
    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)>2);    %最大速度设置，粒子的范围宽度
        v(c) = sign(v(c))*2;   %如果速度大于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(part_size,arr_present(j,1:part_size)); % 适应度（保存至最后一列）

        % 1是适应度评价 2是可行域限制
        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] = min(arr_present(:,end));    % 如果是最小的值为min,相反则为max
    if best<gbest(end) & ( Region_in(arr_present(best_index,:),region) ) % 如果当前最好的结果比以前的好，则更新最优值gbest,如果是最小的值为小于号,相反则为大于号
        gbest = arr_present(best_index,:); % 全局的极值（是第几个粒子？）
    end
    

  
    %------------混沌---------------------------------      
    xlhd = gbest(1:part_size);    
    if(1==2)
       for t=1:20 %次数
        %1生成
              cxl=rand(1,part_size);
              for j=1:part_size
                  if cxl(j)==0
                      cxl(j)=0.1;
                  end
                  if cxl(j)==0.25
                      cxl(j)=0.26;
                  end
                  if cxl(j)==0.5
                      cxl(j)=0.51;
                  end
                  if cxl(j)==0.75
                      cxl(j)=0.76;
                  end
                  if cxl(j)==1
                      cxl(j)=0.9;
                  end
              end
        %2映射
              al=-3;bl=3;
              rxl=al+(bl-al)*cxl;
        %3搜索
              bate = 0.3;
              xlhd=xlhd+bate*rxl;
              if fitness(part_size,xlhd)<gbest(end)
                  gbest(1:part_size)=xlhd;
                  gbest(end)=fitness(part_size,xlhd);
              end
        %4更新
             for j=1:part_size
                  cxl(j)=4*cxl(j)*(1-cxl(j));
             end
       end
    end
    %-------------混沌--------------------------------
    
    
    %当前代的最优粒子的适应度（取自最后一列）保存
    best_record(i) = gbest(end);
end

pso = gbest; % 最优个体
display(gbest);

figure;
plot(best_record);% 最优解与代数的进化关系图


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






%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%返回随机的位置
function ini_present=ini_pos(pop_size,part_size)
ini_present = 5*rand(pop_size,part_size+1);       %初始化当前粒子位置,使其随机的分布在工作空间 
%返回一个随机的矩阵，10*（2+1），最后一列将用来保存适应度

%返回随机的速度
function ini_velocity=ini_v(pop_size,part_size)
ini_velocity =5/2*(rand(pop_size,part_size));     %初始化当前粒子速度,使其随机的分布在速度范围内


%判断是否处于范围内
function flag = Region_in(pos_present,region)
[m n]=size(pos_present); % 1*11  n返回解的维数10
flag=1;
for j=1:n-1
    %pos_present(1:j)是这一行的一个元素
    %flag = flag & ( pos_present(j:1)>=region(j,1) ) & ( pos_present(j:1)<=region(j,2) );
end


%初始化适应度？
function arr_fitness = ini_fit(pos_present,pop_size,part_size)
for k=1:pop_size
    % k行，即每个粒子，的每一列    pos_present(k,1:part_size)表示第k行整体
    % 返回值的第一列保存适应度   arr_fitness 的第一列保存适应度
    arr_fitness(k,1) = fitness(part_size,pos_present(k,1:part_size));  %计算原始种群的适应度    
end

%***************************************************************************
%      计算适应度
%***************************************************************************
function fit = fitness(n,xp)
%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);

%y0=[-10,0]; % 注意：这是基准序列，也就是单个最优的极值
rt1= f1(n,xp);
%rt2= f2(n,xp);
%y1=[rt1,rt2];  % n为粒子维数
%fit=graydegree(2,y0,y1); % 这里关联度在某种意义上就是适应度
fit=rt1;

%目标函数1
function r=f1(n,x)
r=0;
for i=1:n-1
    r=r+(-10)*exp((-0.2)*sqrt(x(i)^2+x(i+1)^2));
end

%目标函数2
function r=f2(n,x)
r=0;
for i=1:n
    r=r+(abs(x(i)))^0.8+5*sin(x(i)^3);
end

%约束函数1
function r=g1(n,x)
r=0;
for i=1:n
    r=0;
end

%约束函数2
function r=g2(n,x)
r=0;
for i=1:n
    r=0;
end

%  灰色关联度计算函数 ( 越大相似性越好 ) 
%  tn目标函数个数   x0基准序列（一组值）   x1贷检（一组值）
function gama = graydegree( tn,y0,y1 )
gama=0;
rou =0.2;
kesa= zeros(tn,1);    
m1= abs(y0(1)-y1(1)) ;
m2= abs(y0(1)-y1(1)) ;
for i=1:tn
    if( abs(y0(i)-y1(i))<m1 )
        m1= abs(y0(i)-y1(i));
    end
    if( abs(y0(i)-y1(i))>m2 )
        m2= abs(y0(i)-y1(i));
    end
end

for i=1:tn
    kesa(i) = ( m1+rou*m2 )/( abs(y0(i)-y1(i)) +rou*m2 );
    gama = gama + kesa(i);
end
gama = gama/tn;


% 可行解的判决函数  gn为约束条件的个数(暂时未用)   n为解（粒子）的维数
function bool = feasible( x,n )
r=0;
%for i=1:gn
    r=max( 0, g1(n,x), g2(n,x) );%判断约束条件
%end
if(r>0)
    bool=0; %不可行解
else
    bool=1; %可行解
end