pso.m

带有测试函数的PSO算法源代码

% File: PSO and Immune 
% Author: xxhalec
% ysu
% m 粒子数
% n 变量维数（暂时为2维，有待改进）
% x m*n 矩阵
% pbx m*n 矩阵 
% v m*n 矩阵
% gbx 1*n 矩阵---全局最优解
% pbx 1*n 矩阵---个体最优解
% pbf m*1 矩阵----个体最优值
% gbf----实数，全局最优解
% k：迭代次数

%格式标准化

clear all;
clc;
format long;

%初始化各个因子
tic
n=2;
m=160;
c1=1.6;
c2=1.6;
w=0.7298;
vmax=100;
xmax=400;
eps=10^(-5); %满足条件的误差
MaxDT=1000; % 最大迭代次数

%产生随机m个粒子,初始化粒子的位置和速度
x=-100+200.*rand(m,n);
v=100.*rand(m,n);

%计算粒子的适应度
for i=1:m
    x(i,3*n+1)=fitness(x(i,:),n);
end
x(:,n+1:2*n)=v;
%按照粒子的适应度升序排列
x(:,3*n+2)=abs(x(:,3*n+1));
B=sortrows(x,3*n+2);


%求个体极值和全局极值
B(:,2*n+1:3*n)=B(:,1:n);
pbx=B(:,2*n+1:3*n);
pbf=B(:,3*n+1);
gbx=B(1,2*n+1:3*n);
gbf=B(1,3*n+1);
xk=B(:,1:n);
v=B(:,(n+1):2*n);

%根据粒子群公式迭代，位置、速度更新
for i=1:m
    v(i,:)=w.*v(i,:)+c1.*rand.*(pbx(i,:)-xk(i,:))+c2.*rand.*(gbx-xk(i,:));
    for j=1:n
        if v(i,j)>vmax
            v(i,j)=vmax;
        elseif v(i,j)<-vmax
            v(i,j)=-vmax;
        end
    end    
    xk(i,:)=xk(i,:)+v(i,:);
     for j=1:n
        if xk(i,j)>xmax
            xk(i,j)=xmax;
        elseif xk(i,j)<-xmax
            xk(i,j)=-xmax;
        end
    end    
end
B(:,1:n)=xk;
B(:,(n+1):2*n)=v;
x=xk;

%开始循环

for k=1:MaxDT
     for i=1:m
            x(i,2*n+1)=fitness(x(i,:),n);
     end
     for i=1:m
         if abs(x(i,2*n+1))<abs(B(i,3*n+1))
             B(i,3*n+1)=x(i,2*n+1);
             B(i,2*n+1:3*n)=x(i,1:n);
         end
     end
     B(:,3*n+2)=abs(B(:,3*n+1));
     B=sortrows(B,3*n+2);
     pbx=B(:,2*n+1:3*n);
     pbf=B(:,3*n+1);
     gbf=B(1,3*n+1);
     gbx=B(1,2*n+1:3*n);
     
     if abs(gbf)<eps
     break
     end
     
     
    
     
     
     xk=B(:,1:n);
     v=B(:,(n+1):2*n);
     for i=1:m
         v(i,:)=w.*v(i,:)+c1.*rand.*(pbx(i,:)-xk(i,:))+c2.*rand.*(gbx-xk(i,:));
         for j=1:n
            if v(i,j)>vmax
                    v(i,j)=vmax;
            elseif v(i,j)<-vmax
                    v(i,j)=-vmax;
            end
         end
         xk(i,:)=xk(i,:)+v(i,:);
         for j=1:n
            if xk(i,j)>xmax
                 xk(i,j)=xmax;
            elseif xk(i,j)<-xmax
                    xk(i,j)=-xmax;
            end
         end    
     end

     B(:,1:n)=xk;
     B(:,(n+1):2*n)=v;
    
     x=xk;   
     
   
    G(k)=gbf;
    G(k)
   
end
%循环结束
G(k)=gbf;


gbx
gbf
k

toc




%画图绘制进化代数和适应度关系曲线图
xx=linspace(1,k,k);
yy=G;
plot(xx,yy,'b-')
hold on
grid on
title('基本粒子群算法进化代数与适应度值关系曲线图')
legend('粒子适应度曲线走势')