pai.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=10;
m=50;
c1=1.6;
c2=1.6;
w=0.7298;

xmax=600;
vmax=100;

eps=10^(-5); %满足条件的误差
MaxDT=1000; % 最大迭代次数

%产生随机m个粒子,初始化粒子的位置和速度
x=-600+1200.*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     
     
     %以下为免疫程序,b---克隆常数，c---变异常数
    b=1;
    c=100;
    E=0;
    for i=1:m
       E=E+B(i,3*n+1);
    end
%开始克隆，克隆个数a由特定公式给出,用ceil取整朝正无穷方向取整,克隆的对象为gbx；
%克隆的对象为x时怎样？？？？需验证？？？
    for i=1:m
        temx=B(i,2*n+1:3*n);
        temf=B(i,3*n+1);
        a=ceil(b*m/i);
        
 %%%克隆出粒子集合
      for j=1:a
          for r=1:n
            cl(j,r)=temx(:,r)+c*temf/E*randn;         
             if cl(j,r)>xmax
                 cl(j,r)=xmax;
             elseif cl(j,r)<-xmax
                    cl(j,r)=-xmax;
             end
          end    
      end
    %选择对抗原具有最大亲和力的克隆抗体
      for j=1:a
        cl(j,n+1)=fitness(cl(j,:),n);
      end
      cl(:,n+2)=abs(cl(:,n+1));
      clB=sortrows(cl,n+2);
    %克隆替换
      if abs(clB(1,n+1))<abs(B(i,3*n+1))
        for r=1:n
            B(i,2*n+r)=clB(1,r);
        end
        B(i,3*n+1)=clB(1,n+1);
      end
      clB=0;
      cl=0;
    end
 %克隆结束                 
     %x(:,(n+1):2*n)=v;
     %x(:,2*n+2)=abs(x(:,2*n+1));
     B(:,3*n+2)=abs(B(:,3*n+1));
     B=sortrows(B,3*n+2);
     
     pbf=B(:,3*n+1);
     gbf=B(1,3*n+1);
     pbx=B(:,2*n+1:3*n);
     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('粒子适应度曲线走势')