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: monkeyLXC72@hotmail.com
//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
