spso.m

许多的PSO优化算法

%测试函数：f(x,y)=100(x^2-y)^2+(1-x)^2,   -2.048<x,y<2.048
%求解函数最小值

global popsize;         %种群规模
%global popnum;         %种群数量
global pop;             %种群
%global c0;             %速度惯性系数,为0—1的随机数
global c1;              %个体最优导向系数
global c2;              %全局最优导向系数
global gbest_x;         %全局最优解x轴坐标
global gbest_y;         %全局最优解y轴坐标
global best_fitness;    %最优解
global best_in_history; %最优解变化轨迹
global x_min;           %x的下限
global x_max;           %x的上限
global y_min;           %y的下限
global y_max;           %y的上限
global gen;             %迭代次数
global exetime;         %当前迭代次数
global max_velocity;    %最大速度

%初始化**********************************************************************************************************************************

gen=100;                      %设置进化代数
popsize=30;                   %设置种群规模大小
best_in_history(gen)=inf;     %初始化全局历史最优解
best_in_history(:)=inf;       %初始化全局历史最优解
max_velocity=0.3;             %最大速度限制
best_fitness=inf;
%popnum=1;                    %设置种群数量

pop(popsize,8)=0;   
%初始化种群,创建popsize行5列的0矩阵
%第1列为x轴坐标，第2列为y轴坐标
%第3列为x轴速度分量，第4列为y轴速度分量
%第5列为个体最优位置的x轴坐标，第6列为个体最优位置的y轴坐标
%第7列为个体最优适值，第8列为当前个体适应值

for i=1:popsize
    pop(i,1)=4*rand()-2;            %初始化种群中的粒子位置，值为-2—2，步长为其速度
    pop(i,2)=4*rand()-2;            %初始化种群中的粒子位置，值为-2—2，步长为其速度
    pop(i,5)=pop(i,1);              %初始状态下个体最优值等于初始位置
    pop(i,6)=pop(i,2);              %初始状态下个体最优值等于初始位置
    pop(i,3)=rand()*0.02-0.01;      %初始化种群微粒速度，值为-0.01—0.01，间隔为0.0001
    pop(i,4)=rand()*0.02-0.01;      %初始化种群微粒速度，值为-0.01—0.01，间隔为0.0001
    pop(i,7)=inf;
    pop(i,8)=inf;
end

c1=2;
c2=2;
x_min=-2;
y_min=-2;
x_max=2;
y_max=2;

gbest_x=pop(1,1);   %全局最优初始值为种群第一个粒子的位置
gbest_y=pop(1,2);

%初始化*********************************************************************************************************************************

for exetime=1:gen  
    
%实时输出结果***************************************************************************************************************************

%输出当前种群中粒子位置
subplot(1,2,1);
for i=1:popsize
    plot(pop(i,1),pop(i,2),'b*');
    hold on;
end

plot(gbest_x,gbest_y,'r.','markersize',20);axis([-2,2,-2,2]);
hold off;

subplot(1,2,2);
axis([0,gen,-0.00005,0.00005]);

if exetime-1>0
    line([exetime-1,exetime],[best_in_history(exetime-1),best_fitness]);hold on;
end 
%实时输出结果****************************************************************************************************************************

%适值计算*******************************************************************************************************************************

% 测试函数为f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048  计算适应值并赋值
for i=1:popsize
    pop(i,8)=100*(pop(i,1)^2-pop(i,2))^2+(1-pop(i,1))^2;
    if pop(i,7)>pop(i,8)                    %若当前适应值优于个体最优值，则进行个体最优信息的更新
        pop(i,7)=pop(i,8);                  %适值更新
        pop(i,5:6)=pop(i,1:2);              %位置坐标更新
    end
end

%计算完适应值后寻找当前全局最优位置并记录其坐标
if best_fitness>min(pop(:,7))
    best_fitness=min(pop(:,7));                      %全局最优值
    gbest_x=pop(find(pop(:,7)==min(pop(:,7))),1);    %全局最优粒子的位置   
    gbest_y=pop(find(pop(:,7)==min(pop(:,7))),2);
end

best_in_history(exetime)=best_fitness;               %记录当前全局最优
    
%适值计算*******************************************************************************************************************************

%粒子群速度与位置更新********************************************************************************************************************

%更新粒子速度
for i=1:popsize
    pop(i,3)=rand()*pop(i,3)+c1*rand()*(pop(i,5)-pop(i,1))+c2*rand()*(gbest_x-pop(i,1));    %更新速度
    pop(i,4)=rand()*pop(i,4)+c1*rand()*(pop(i,6)-pop(i,2))+c2*rand()*(gbest_x-pop(i,2)); 
    if abs(pop(i,3))>max_velocity
        if pop(i,3)>0
            pop(i,3)=max_velocity;
        else
            pop(i,3)=-max_velocity;
        end
    end
    if abs(pop(i,4))>max_velocity
        if pop(i,4)>0
            pop(i,4)=max_velocity;
        else
            pop(i,4)=-max_velocity;
        end
    end
end

%更新粒子位置
for i=1:popsize
    pop(i,1)=pop(i,1)+pop(i,3);
    pop(i,2)=pop(i,2)+pop(i,4);
end

pause(0.01);
end

clear i;
clear exetime;
clear x_max;
clear x_min;
clear y_min;
clear y_max;