[求助] 期刊文章《基于空间缩放和吸引子的粒子群优化算法》的算法实现问题

clear all;
clc;
format short;

D=30;
P=15;
iter_max=5000;
Xmin=-100;
Xmax=100;
func_num=2;     
pop=rand(P,D)*200-100;
pop_val=ones(1,P);
global initial_flag
initial_flag=0;

SzAPSO(pop,pop_val,D,P,iter_max,Xmin,Xmax,func_num);

function [Pg_val,Pg,fitcount] = SzAPSO(pop,  pop_val, D, P,     iter_max, Xmin,  Xmax,  func_num)
                                       %粒子 适应值   维 粒子数 最大迭代  最小值 最大值 函数代号

popbest = pop;  %记录局部最优解
Pg = pop(1,:);
Pg_val = benchmark_func(pop(1,:), func_num);

for i = 1:P                         %初始化粒子最优和全局最优
    pop_val(i) = benchmark_func(pop(i,:), func_num);
    if pop_val(i) < Pg_val
        Pg_val = pop_val(i);
        Pg = pop(i,:);
        end
end

fitcount = 0;
k = 100;

while fitcount<iter_max
    z = k * rand;
    for i = 1:P
        if abs(pop_val(i)-Pg_val) < 1e-6               %判断粒子是不是最优，计算吸引子Pa
            Pa = Pg * (rand+0.5) + (rand*2-1);          %围绕全局最优抖动
        else
            Pa = popbest(i)+(Pg-popbest(i)) * rand;     %局部最优和全局最优加权平均
        end
        
        for j = 1:D
            if Pa(j) > Xmax
                Pa(j) = Xmax;
                if pop(i,j) > Pa(j)
                    pop(i,j) = Xmax;
                else
                    pop(i,j) = update(pop(i,j),Pa(j),z);
                end
            elseif Pa(j) < Xmin
                Pa(j) = Xmin;
                if pop(i,j) < Pa(j)
                    pop(i,j) = Xmin;
                else
                    pop(i,j) = update(pop(i,j),Pa(j),z);
                end
            else
                pop(i,j) = update(pop(i,j),Pa(j),z);
            end
        end
        
        pop_val(i) = benchmark_func(pop(i,:), func_num);        %更新全局最优
        
        if pop_val(i) < benchmark_func(popbest(i), func_num)  %更新局部最优
            popbest(i) = pop(i);
        end
        
        if pop_val(i) < Pg_val      %更新全局最优
            Pg_val = pop_val(i);
            Pg = pop(i,:);
        end
        
    end

    fprintf('%f ,  %d\n\n',Pg_val,fitcount);
    if Pg_val<1e-6      %满足精度结束
        break;
    end
    fitcount = fitcount + 1;
%       if fitcount == 1000
%          keyboard;
%       end
end

disp('************************function*****************************');
disp(func_num);
disp('The best result is:')
disp(Pg_val);
disp(fitcount);
disp('*************************************************************');
keyboard;
clear all;
end

function num = update(Xid,Pa,z)     %更新后会有超界的情况(但是取模后在加上Pa应该是一定在取值空间内才对)

if Xid < Pa
    num = mod((Xid-Pa),((Pa - -100)*z))/z+Pa;
else
    num = mod((Xid-Pa),((Pa * -1 +100)*z))/z+Pa;
end

end

function num = schwefel_PSO(pop,b)

num = 0;
a = length(pop);
temp(a) = 0;
for i = 1:a
    for j = 1:i
        temp(i) = temp(i) + pop(j);
    end
    num = num + temp(i)^2;
end

end