liziqun.txt

基于Scilab的粒子群算法代码

function [f_opt,x_opt]=PSO(fun,dim,itval_min,itval_max,num,typ,time,v_step,c0_strth,c1_strth,c2_strth,figon)
//***********************************************
//PSO algorithm for finding the maximal/minimal value of a function defined on n dimention space
//input parameter
//           fun: <function>  optimization function which should be in the form of [y]=fun(x)
//                            x=[x(1) x(2) x(3) ... ect.] where x(n) is the nth variation
//           dim: <integer>   dimention of the assemble of definition of the function
//     itval_min:  <array>    the minimum boundaries of each dimention
//     itval_max:  <array>    the maximum boundaries of each dimention
//         [num]: <integer>   the number of particals
//         [typ]:    <1>      find the maximal value of the function
//                 <other>    find the minimal value of the function
//      [v_step]:  <float>    v_step times defaut reference velociy depend on the boundary of each dimention
//                 <array>    reference velocity for each dimention
//    [c0_strth]:  <[0,1]>    the strength parameter for the present velocity
//    [c1_strth]:  <[0,1]>    the strength parameter for the local attractors
//    [c2_strth]:  <[0,1]>    the strength parameter for the global attractor
//       [figon]:    <1>      show figure 
//output parameter
//         x_opt:  <array>    the coordinates of the max/min value of fun
//         f_opt:  <float>    the value of fun(x)
//
//exemple
//  deff('[y]=fun(x)','y=100*(x(:,1)-x(:,2).^2)^2+(x(:,2)-1)^2');
//  dim=2;
//  ivmax=[5,5];
//  ivmin=[-5,-5];
//  [x]=PSO(fun,dim,ivmin,ivmax,20,0,500,1,0.99,0.8,0.8,1)
//  [x]=PSO(fun,dim,ivmin,ivmax,20,0,500) //for fast
//
//code by    
//Liu. Xinchang
//Ecole Centrale de Pekin (BeiHang University)
//China
//email: [e-mail]
//reference: the cours of modeling by Hu Baogang
//25 Mar, 2008.
//***********************************************

//--  Check  ----------------------------------
f_opt=0;
x_opt=0;
if size(itval_max,2)<>dim|size(itval_min,2)<>dim
  if size(itval_max,1)<>dim|size(itval_min,1)<>dim    
  printf("Err: interval size error!")
    return;
  else
    itval_max=itval_max';
    itval_min=itval_min';
  end
end
if min(itval_max-itval_min)<0
    printf("Err: interval max-min error!");
  return;
end

//--  initialization  -------------------------
if ~exists('typ')
  printf("PSO_output: Find the maximum value of the function\n");
elseif typ==1
  printf("PSO_output: Find the maximum value of the function\n");
else
  typ=-1
  printf("PSO_output: Find the minimum value of the function\n");
end
if ~exists('num')
  num=dim^2+1;
  printf("PSO_output: Use defaut set of num\n");
elseif num<=1
  num=dim^2+1;
  printf("Warning: num is too small\n");
else
  num=sum(num);
end
if ~exists('time')
  time=1000;
  printf("PSO_output: Use defaut set of time\n");
end
if ~exists('v_step')
  v_step=(itval_max-itval_min)/1000;
  printf("PSO_output: Use defaut set of v_step\n");
elseif size(v_step)==[1,1] 
  v_step=v_step.*(itval_max-itval_min)/1000;
elseif size(v_step)==[dim,1]
  v_step=v_step';
elseif size(v_step)<>[1,dim]
  v_step=ones(1,dim)*v_step(1);
  printf("Warn: length of v_step is not suitable\n");
end
if ~exists('c0_strth')
  c0_strth=0.99;
  printf("PSO_output: Use defaut set of c0_strth\n");
end
if ~exists('c1_strth')
  c1_strth=mean(v_step.^2)^0.5*0.8;
  printf("PSO_output: Use defaut set of c1_strth\n");
else
  c1_strth=c1_strth*mean(v_step.^2)^0.5;
end
if ~exists('c2_strth')
  c2_strth=mean(v_step.^2)^0.5*0.8;
  printf("PSO_output: Use defaut set of c2_strth\n");
else
  c2_strth=c2_strth*mean(v_step.^2)^0.5;
end
if ~exists('figon')
  figon=0
end  
position=rand(num,dim)*diag(itval_max-itval_min)+ones(num,dim)*diag(itval_min);
velocity=(rand(num,dim)-0.5)*diag(v_step);
if time<0|time>1000000
  time=1000;
end
t=0;
for j=1:num
  opt_local(j,:)=typ*fun(position(j,:));
end
pst_opt_loc=position;

//- algorithme  -------------------------------
while t<time do
  t=t+1;
  for j=1:num
    [opt_local(j,:),i(j)]=max(typ*fun(position(j,:)),opt_local(j));
  end
  pst_opt_loc=diag(2-i)*ones(num,dim).*position+diag(i-1)*ones(num,dim).*pst_opt_loc;
  [opt_global,i]=max(opt_local);
  pst_opt_glo=pst_opt_loc(i,:);
  v_loc=pst_opt_loc-position;
  v_glo=ones(num,dim)*diag(pst_opt_glo)-position;
  velocity=c0_strth*velocity+c1_strth*(diag((sum(v_loc.^2,2)+0.0001).^(-0.5))*v_loc).*rand(num,dim)+c2_strth*(diag((sum(v_glo.^2,2)+0.0001).^(-0.5))*v_glo).*rand(num,dim)
  position=position+velocity;
  position=max(position,ones(num,dim)*diag(itval_min)); //check boundaries
  position=min(position,ones(num,dim)*diag(itval_max));
  if figon==1&size(position,2)>1
    plot(position(:,1),position(:,2),'r.','marksize',1)
  end
end

//--  output  ---------------------------------
x_opt=pst_opt_glo;
f_opt=fun(x_opt);
printf("PSO_output: the result is (%s), the value is %f\n",strcat(string(x_opt),','),f_opt);
return;
endfunction