qq11eeee.m

this modified PSO algorithm can solve for constrained optimization problem

 function [pso F] = pso_2D_new()

clc;
clear;close all;
N = 10;          %   N 种群大小
D = 38;          %   D 粒子维数,                                                                      ** =n-D
r=5;             %   移动锚节点的发射半径，单位为米
L=20;            %   监测区域宽与高，单位为米
S=30;            %   未知节点个数
M = [1 3;4 7;3 10;7 2;9 8;3 16;14 17;16 13;20 10;17 18;
     5 9;6 6;3 5;8 17;4 20;5 18;7 14;10 12;11 14;13 15;
     15 14;16 20;13 8;14 11;12 5;15 3;18 7;16 9;11 10;8 11];
max_gen = 100;          %   max_gen 最大迭代次数
%region=zeros(D,2);  %   设定搜索空间范围
region=[0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20;0,20;0,20;
    0,20;0,20];          %                                                                                      **每一维设定不同范围


rand('state',sum(100*clock));   %   重置随机数发生器状态
pop = ini_pos(N,D);   %   present 当前位置,随机初始化,rand()的范围为0~1

v=ini_v(N,D);             %   初始化当前速度

gbest = zeros(1,D);            %   gbest 当前搜索到的最小的值
pbest = zeros(N,D);      %   pbest 粒子以前搜索到的最优值，最后一列包括这些值的适应度
w_max = 0.9;                            %   w_max 权系数最大值
w_min = 0.4;
v_max = 6;             %                                                                                          **最大速度,为粒子的范围宽度
c1 = 2;                   %   学习因子
c2 = 2;                   %   学习因子
best_record = zeros(1,max_gen);     %   best_record记录发生max_gen次迭代得到的最好的粒子的适应度。
%  ————————————————————————
%   计算原始种群的适应度,及初始化
%  ————————————————————————
for j=1:N
    f_pop(j)=fitness(pop(j,:),D,r,L);     %将N个粒子的适应值赋给f_pop
end
%pbest = pop;                                        %初始化各个粒子最优值

f_pbest=f_pop;                %f_pbest各个粒子的个体极值(最好适应值)
[best_value best_index] = min(f_pbest);         %初始化全局最优，即适应度为全局最小的值，根据需要也可以选取为最大值
gbest = pop(best_index,:);    %最小适应值所在行的所有列赋给gbest，
f_gbest=f_pbest(best_index);  %f_gbest为全局最优适应值

x=[0:0.01:20];
y=[0:0.01:20];
%z=@(x,y) (1-x).^2 - (y+1).^2+(1+x).^2-(y+2).^2 ;
%z=@(x,y) (3+x).^2 - (y+1).^2 ;
z=@(x,y) sum(x-y)^2;
h=waitbar(0,'Please wait...');
for i=1:max_gen
    grid on;
    ezmesh(z),hold on,grid on,plot3(pop(:,1),pop(:,2),pop(:,3),pop(:,4),pop(:,5),pop(:,6),pop(:,7),pop(:,8),pop(:,9),pop(:,10),pop(:,11),pop(:,12), pop(:,13),pop(:,14),pop(:,15),pop(:,16),pop(:,17),pop(:,18),pop(:,19),pop(:,20),pop(:,21),pop(:,22),pop(:,23),pop(:,24),pop(:,25),pop(:,26),pop(:,27),pop(:,28),pop(:,29),pop(:,30),pop(:,31),pop(:,32),pop(:,33),pop(:,34),pop(:,35),pop(:,36),pop(:,37),pop(:,38),f_pop','r*'),hold off;
    %ezmesh(z),hold on,grid on,plot3(pop(:,1),pop(:,2),f_pop','r*'),hold off;
    F(i)=getframe;
    
    %   确定是否对打散已经收敛的粒子群——————————————————————————————
    reset = 0;          %   reset = 1时设置为粒子群过分收敛时将其打散，如果＝1则不打散
    if reset==1
        bit = 1;
        for k=1:D
            bit = bit&(range(pop(:,k))<0.1);
        end
        if bit==1       %   bit=1时对粒子位置及速度进行随机重置
            pop = ini_pos(N,D);   %   present 当前位置,随机初始化
            v = ini_v(N,D);           %   速度初始化
            for j=1:N
                f_pop(j)=fitness(pop(j,:),D,r,L);
            end
            warning('粒子过分集中！重新初始化……');      %   给出信息
            display(i);
        end
    end

    w = w_max-(w_max-w_min)*i/max_gen;
    %根据基本粒子群迭代公式进行一次进化迭代产生中间子代pop(j,1:D)
    for j=1:N
        v(j,:) = w.*v(j,:)+c1.*rand.*(pbest(j,1:D)-pop(j,1:D))...
            +c2.*rand.*(gbest(1:D)-pop(j,1:D));                        %  粒子速度更新 (a)

        %   判断v的大小，限制v的绝对值小于5————————————————————————————
       % c = find(abs(v)>v_max);                                                                                              %**最大速度设置，粒子的范围宽度
       % v(c) = sign(v(c))*v_max;   %如果速度大于6则，速度为6

        pop(j,1:D) = pop(j,1:D)+v(j,1:D);              %  粒子位置更新 (b)
        f_pop(j) = fitness(pop(j,1:D),D,r,L);
        if violent(M,pop(j,1:D),S,D,r)>=3
           if (f_pop(j)<f_pbest(j))&(Region_in(pop(j,:),region))     %   根据条件更新pbest,如果是最小的值为小于号,相反则为大于号
             pbest(j,:) = pop(j,:);
             f_pbest(j)=f_pop(j);
           end
        end

    end

    [best_value best_index] = min(f_pbest);                                      %   如果是最小的值为min,相反则为max
    
    if (best_value<f_gbest)&(Region_in(pop(best_index,:),region))                  %   如果当前最好的结果比以前的好，则更新最优值gbest,如果是最小的值为小于号,相反则为大于号
        gbest = pop(best_index,:);
        f_gbest=f_pbest(best_index);
    end

    best_record(i)=f_gbest;
    waitbar(i/max_gen);
end
close(h);
pso = gbest;
display(f_gbest);
figure;plot(gbest(1:2:end-1),gbest(2:2:end),'go'); 


%名称：粒子违反约束程度函数
%输出：粒子违反约束条件的程度（若违反输出大于，不违反输出等于零）
%形式：保存为M文件被调用

%约束函数
function  number=violent(M,present,S,D,r) 
for k=1:S
    number=0;
    for t=1:2:D-1
      if(M(k,1)-present(t))^2+(M(k,2)-present(t+1))^2<=r^2
          number=number+1;
      end
    end
end  


%   ***************************************************************************
%      计算适应度
%   ***************************************************************************
function fit = fitness(present,D,r,L)
%fit=(3+present(:,1)).^2 - (present(:,2)+1).^2 +(3+present(:,3)).^2-(present(:,4)+1).^2
   % +(3+present(:,5)).^2 - (present(:,6)+1).^2 +(3+present(:,7)).^2-(present(:,8)+1).^2
   % +(3+present(:,9)).^2 - (present(:,10)+1).^2 +(3+present(:,11)).^2-(present(:,12)+1).^2
   % +(3+present(:,13)).^2 - (present(:,14)+1).^2 +(3+present(:,15)).^2-(present(:,16)+1).^2
   % +(3+present(:,17)).^2 - (present(:,18)+1).^2 +(3+present(:,19)).^2-(present(:,20)+1).^2
   % +(3+present(:,21)).^2 - (present(:,22)+1).^2 +(3+present(:,23)).^2-(present(:,24)+1).^2
    %+(3+present(:,25)).^2 - (present(:,26)+1).^2 +(3+present(:,27)).^2-(present(:,28)+1).^2
    %+(3+present(:,29)).^2 - (present(:,30)+1).^2 +(3+present(:,31)).^2-(present(:,32)+1).^2
    %+(3+present(:,33)).^2 - (present(:,34)+1).^2 +(3+present(:,35)).^2-(present(:,36)+1).^2
    %+(3+present(:,37)).^2 - (present(:,38)+1).^2 ;
    
    
%适应函数    
for j=1:D/2-1
    count=0;
    fit=0;
    average=pi*r*r*D/2*L*L;
      for i=j+1:D/2
        if (present(j)-present(2*i-1))^2+(present(j+1)-present(2*i))^2<=r^2
           count=count+1;
        end
      end
      fit=fit+(count-average)^2;
end



function ini_position=ini_pos(N,D)
ini_position = 20*rand(N,D);        %初始化当前粒子位置,使其随机的分布在工作空间                         %** 6即为自变量范围


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


function flag=Region_in(present,region)
flag=1;
len=length(present);
for k=1:len
    flag=flag&(present(k)>=region(k,1))&(present(k)<=region(k,2));
end