parapso.m

智能微粒群为最小二乘支持向量机调参的示例程序（LSSVM+PSO）

function C = parapso(trnX, trnY, tstX, tstY, PSOoptions)
%SIMPLEGENETICALGORITHM      To optimize parameters by using the PSO algorithm.
ker='rbf'; 
global p1; 
nps = PSOoptions(1);  % number of particles
iter = PSOoptions(2); % number of algorithm iterations
Wmin = PSOoptions(3); % minimum inertial weight
Wmax = PSOoptions(4); % maximum inertial weight
VCmin = PSOoptions(5); % minimum change in model parameter
VCmax = PSOoptions(6); % maximum change in model parameter
Vpmin = PSOoptions(7); % minimum change in model parameter
Vpmax = PSOoptions(8); % maximum change in model parameter
c1 = PSOoptions(9);   % parameter in PSO methodology
c2 = PSOoptions(10);   % parameter in PSO methodology
Cmin = PSOoptions(11);% minimum model parameter value
Cmax = PSOoptions(12);% maximum model parameter value
p1min = PSOoptions(13);% minimum model parameter value
p1max = PSOoptions(14);% maximum model parameter value
%initialize the particles and their velocity components
% for i=1:5
%     S(i) = zeros(20 ,1);
%     SY(i) = zeros(20 ,1); 
%     trnX(i) = zeros(80 ,1); 
%     trnY(i) = zeros(80 ,1);
% end
% S(1)=X(1:20,1);SY(1)=Y(1:20,1);
% trnX(1)=X(21:100,1);trnY(1)=Y(21:100,1);
% S(2)=X(21:40,1);SY(2)=Y(21:40,1);
% trnX(2)=X(1:20,1)+X(41:100,1);trnY(2)=Y(1:20,1)+Y(41:100,1);
% S(3)=X(41:60,1);SY(3)=Y(41:60,1);
% trnX(3)=X(1:40,1)+X(61:100,1);trnY(3)=Y(1:40,1)+Y(61:100,1);
% S(4)=X(61:80,1);SY(4)=Y(61:80,1);
% trnX(4)=X(1:60,1)+X(81:100,1);trnY(4)=Y(1:60,1)+Y(81:100,1);
% S(5)=X(81:100,1);SY(5)=Y(81:100,1);
% trnX(5)=X(1:80,1);trnY(5)=Y(1:80,1);

   %t == 1;
    for i = 1:nps
     Ctemp(i) = Cmax*rand(1);
     ptemp(i) = p1max*rand(1);
     p1 = ptemp(i);
     VC(i) = (VCmax-VCmin)*rand(1) + VCmin;
     Vp(i) = (Vpmax-Vpmin)*rand(1) + Vpmin;
     pbestC(i) = Ctemp(i); 
     pbestp1(i)= ptemp(i);
     temp = Ctemp(i);
     % for j=1:5
        [alpha,bias]=calr(trnX,trnY,ker,temp);%最小二乘LU
        funY=output(trnX,tstX,ker,alpha,bias);
        RMSE=norm(funY-tstY)/sqrt(length(tstY)); 
     % end 
     pbestfit(i)=RMSE;
    end
    gbestfit = min(pbestfit);
    for i=1:nps
        if pbestfit(i)==gbestfit
            gbestC = Ctemp(i);
            gbestp1 = ptemp(i);
        end
    end

 
  for k = 2:iter
    for i = 1:nps
     W = Wmax - (k-1)*(Wmax-Wmin)/iter; % compute inertial weight  
     VC(i) = W* VC(i)+c1*rand(1)*(pbestC(i)-Ctemp(i))+c2*rand(1)*(gbestC-Ctemp(i));
     Vp(i) = W* Vp(i)+c1*rand(1)*(pbestp1(i)-ptemp(i))+c2*rand(1)*(gbestp1-ptemp(i));
      if VC(i) > VCmax, Ctemp(i) = Ctemp(i) + VCmax; end
      if VC(i) < VCmin, Ctemp(i) = Ctemp(i) + VCmin; end
      if VC(i) > VCmin & VC(i) < VCmax, Ctemp(i) = Ctemp(i) + VC(i); end
      if Vp(i) > Vpmax, ptemp(i) = ptemp(i) + Vpmax;end
      if Vp(i) < Vpmin, ptemp(i) = ptemp(i) + Vpmin;end
      if Vp(i) > Vpmin & Vp(i) < Vpmax, ptemp(i) = ptemp(i) + Vp(i);end
      if Ctemp(i)< Cmin,Ctemp(i)= Cmin;end
      if Ctemp(i)> Cmax,Ctemp(i)= Cmax;end
      if ptemp(i)< p1min,ptemp(i)= p1min;end
      if ptemp(i)> p1max,ptemp(i)= p1max;end
      temp = Ctemp(i);
      p1 = ptemp(i);
      % for j=1:5
        [alpha,bias]=calr(trnX,trnY,ker,temp);%最小二乘LU
        funY=output(trnX,tstX,ker,alpha,bias);
        RMSE=norm(funY-tstY)/sqrt(length(tstY)); 
      % end 
      if pbestfit(i)> RMSE,
          pbestfit(i)= RMSE;
          pbestC(i)=Ctemp(i);
          pbestp1(i)=ptemp(i);
      end
      if gbestfit >  pbestfit(i),
          gbestfit =  pbestfit(i);
          gbestC = Ctemp(i);
          gbestp1 = ptemp(i);
      end
     end
   end
      C = gbestC;
      p1 = gbestp1;