带收缩因子的PSO优化算法

%--ARPOS
%--带有收缩因子，并保证种群的多样性的微粒群算法解决有约束的情况--

 

%清除屏幕
clc
clear

%-----------------------------参数设置-------------------------------------
c1 = 2;              % 学习因子1 ，一般在[0,2]
c2 = 2;              % 学习因子2 ，一般在[0,2]
% c1 = 2.04344;              %学习因子1 ，一般在[0,2]
% c2 = 0.94874;              %学习因子2 ，一般在[0,2]
k1  = 0.7298;             % 收缩因子 
Dimension = 2;            % 搜索空间维数（未知数个数） 
Popsize = 20;             % 初始化群体个体数目 
MaxDT = 100;              % 最大迭代次数 
DivH      =   0.25;       % 最大多样性系数
DivL      =   0.0005;   % 最小多样性系数

itera = 1;                % 迭代次数

%--------------初始化种群的个体(可以在这里限定位置和速度的范围)--------------- 

% 1）随机初始化位置
x1 = rand(Popsize,1)*(100-13) + 13;         % x1 = [13,100]
x2 = rand(Popsize,1)*100;                   % x2 = [0, 100]

x1minmax = [13 100];
x2minmax = [0 100];
sumL = sqrt((x1minmax(2)-x1minmax(1))^2 + (x2minmax(2)-x2minmax(1))^2);

% 2）随机初始化速度
v1max = 0.5*( max(x1) - min(x1));
v2max = 0.5*( max(x2) - min(x2));
v1 = rand(Popsize,1) * v1max;
v2 = rand(Popsize,1) * v2max;

fit = fitness(x1,x2,itera,Popsize);    % 计算各微粒的适应度
Pbest = [x1 x2];           % 初始化个体（微粒）最佳位置
Fbest = fit;               % 初始化个体（微粒）最佳位置时的适应度
[minfit, indfit] = min(Fbest);   % 寻找群体全局适应度最佳的个体（微粒）[适应度，微粒号]
PGbest = [x1(indfit) x2(indfit)];    % 全局适应度最佳的个体（微粒）位置
FGbest = minfit;                    % 全局适应度最佳的个体（微粒）适应度

dir = 1;
%--------------进入主要循环，按照公式依次迭代，直到满足精度要求--------------- 
for itera=2:MaxDT
    
     div1 = diversity(Popsize, Dimension, Pbest, sumL);
    if dir > 0 && div1 < DivL
        dir = -1; 
    elseif dir < 0  && div1 > DivH
        dir = 1;
    end
    
    dir1(itera) = dir;

    w = 1.2-1.1*itera/MaxDT;% 惯性权重，一般取[0,1.4]，但[0.8,1.2]收敛速度更快 
%-------------更新粒子飞行速度-------------
v1 = k1*(w * v1 + dir*(c1 * rand * (Pbest(:,1) - x1) + c2 * rand * (PGbest(1) - x1)));    % 更新 v1
v2 = k1*(w * v2 + dir*(c1 * rand * (Pbest(:,2) - x2) + c2 * rand * (PGbest(2) - x2)));    % 更新 v2
    
%-------------限制粒子飞行速度-------------    
    v1(find( v1 > v1max)) = v1max;  v1(find( v1 < -v1max)) = -v1max;
    v2(find( v2 > v2max)) = v2max;  v2(find( v2 < -v2max)) = -v2max;
    
%-------------更新粒子位置 -------------       
    x1 = x1 + v1;            % 更新 x1
    x2 = x2 + v2;            % 更新 x2
    
%-------------限制粒子位置-------------       
    x1(find( x1 > 100 ))  = 100; x1(find( x1 < 13 )) = 13;
    x2(find( x2 > 100 ))  = 100; x2(find( x2 < 0 )) =  0;
        
    [new_fit]= fitness(x1, x2, itera, Popsize);     % 计算微粒的适应度
    
%-------------更新个体（微粒）最佳适应度和最佳位置-------------      
    for i = 1 : Popsize
        if  new_fit(i) < Fbest(i)
            Pbest(i,:) = [x1(i) x2(i)];
            Fbest(i) = new_fit(i);
        end
    end

    [minfit(itera), indfit] = min(new_fit);   % 在本次迭代中，寻找最佳适应度值和微粒号

%-------------更新全局最佳适应度和最佳位置-------------    
    if minfit(itera) < PGbest
        PGbest = [x1(indfit) x2(indfit)];    % updating gworst
        FGbest = minfit(itera);             % updating gworst
    end;
    
%-------------最佳适应度和最佳位置-------------      
    PG_BESTx(itera,:) = [x1(indfit) x2(indfit)];
    FG_BEST(itera) =  opti(PGbest(1),PGbest(2));
end    
 PGbest
%  FGbest
% PG_BESTx
min(FG_BEST)
%-------------绘图-------------    
figure(1)
grid on
t = 1:MaxDT;

subplot(2,2,1);
plot(t,FG_BEST,'b');
grid on
title('函数最优值与迭代次数的关系')
xlabel('迭代次数 i');
ylabel('函数最优值 FG-BEST');

subplot(2,2,2);
plot(PG_BESTx(t,1),PG_BESTx(t,2),'*b');
hold on
plot(PG_BESTx(MaxDT,1),PG_BESTx(MaxDT,2),'dr','LineWidth',2);
grid on
title('最佳微粒位置')
xlabel('微粒位置x(1)');
ylabel('微粒位置x(2)');

subplot(2,2,3);
plot(t,dir1,'db');
grid on
title('diversity')
xlabel('迭代次数 i');
ylabel('dir值');    
D132