/*Before running this script, you should run the script "Test functions" to initialize the functions we'll be
working on*/
/*-------------------------------------------------------PSO code----------------------------------------------------*/
/* PSO parameters */
%let nbx = 2; /* Number of variables */
%let p = 20; /* Number of particles */
%let g = 19 ;/* Number of relations between particles : we take p-1 for the global topology */
%let t = 1000; /* Number of iterations */
%let UB = 5; /* Upper bound on solution */
%let LB = -5; /* Lower bound on solution */
%let vmax =(&UB-&LB)/2; /* Maximum velocity */
%let c1 = 2.05; /* Weight on pull toward best particle position */
%let c2 = 2.05; /* Weight on pull toward best global position */
%let func=rosenbrock; /*The function to minimize*/
%let r1=uniform(0); /* Random variable */
%let r2=uniform(0); /* Random variable */
%let c = &c1 + &c2;/* Acceleration coefficient */
%let chi =2/((&c-2)+sqrt(&c**2-4*&c)); /* Constriction coefficient */
%let w=0.5*uniform(0)+0.5; /* Inertia weight */
/* PSO code */
data solution_canonical;
/* Define arrays for particle position(=x), velocity(=v), bests(pbest and gbest)
and one to store the pbests of the objective function */
array x[&nbx,&p];
array v[&nbx,&p];
array pbest[&nbx,&p];
array gbest[&nbx];
array fpbest[&p];
array n[&nbx,&g];
array nbest[&nbx,&g];
array fnbest[&g];
/* Initialize arrays with random particle locations and velocities */
do i = 1 to &p;
do z = 1 to &g;
do j = 1 to &nbx;
x[j,i] = &LB + (&UB - &LB)*&r1; /* Random particles locations */
n[j,z] = &LB + (&UB - &LB)*&r1; /* Random neighbors locations */
v[j,i] = -&vmax +(2*&vmax)*&r1; /* Random particles velocities */
pbest[j,i] = x[j,i];
nbest[j,z]= n[j,z];
fpbest[i] = %f(%str(pbest[:,i]), func=&func);
fnbest[z]= %f(%str(nbest[:,z]), func=&func);
end;
if i = 1 and z = 1 then do;
do j = 1 to &nbx;
gbest[j] = x[j,i];
end;
/* Store the value of the objective function at the global best */
fgnbest = fnbest[z];
fgbest = fpbest[i];
end;
else do;
if fnbest[z] < fgnbest then do;
do j = 1 to &nbx;
gbest[j] = n[j,z];
end;
if fpbest[i] < fnbest[z] then do;
do j = 1 to &nbx;
gbest[j] = x[j,i];
end;
if fpbest[i] < fgbest then do;
do j = 1 to &nbx;
gbest[j] = x[j,i];
end;
end;
fgbest = fpbest[i];
end;
end;
end;
end;
/* Here is the fun part : Update particle locations and velocities to let them wander*/
/* For each iteration, particle and variable */
if t=&t then stop;
do i = 1 to &p;
do z=1 to &g;
do j = 1 to &nbx;
do until ( x[j,i]-x[j,i-1] = 1e-7);
/* Compute new velocity and restrict it within minima and maxima */
v[j,i] = &chi*(&w*v[j,i] + &c1 * &r1 * (pbest[j,i] - x[j,i]) + &c2 * &r2 * (gbest[j] - x[j,i]));
v[j,i] = min(&vmax, max(-&vmax, v[j,i]));
/* Move particle to new position and restrict them to the UP and LB */
x[j,i] = x[j,i]+ v[j,i];
x[j,i] = min(&UB, max(&LB, x[j,i]));
end;
end;
/* Calculate the best particle in each neighborhood and THEN update particle and global best positions (if new best found) */
fx= %f(%str(x[:,i]), func=&func);
if fx < fnbest[z] then do;
do j = 1 to &nbx;
nbest[j,z] = n[j,z];
fpbest[i] = fnbest[z];
end;
if fx < fpbest[i] then do;
do j = 1 to &nbx;
pbest[j,i] = x[j,i];
fpbest[i] = fx;
end;
if fx < fgbest then do;
do j = 1 to &nbx;
gbest[j] = x[j,i];
end;
fgbest = fx;
end;
end;
*end;
end;
end;
end;
end;
/* Record the iterations and global best position as the output solution */
/*do j = 1 to &nbx;
do i = 1 to &p;
var_iteration_indice = "x" || strip(i) || strip(j) ;
iteration_solution = pbest[j,i];
output;
end;
end;
keep var_iteration_indice iteration_solution;*/
do j = 1 to &nbx;
var_indice = "x" || strip(j) ;
global_solution = gbest[j];
output solution_canonical;
end;
keep var_indice global_solution ;
run;
proc transpose data=Solution_canonical out=newt;
var global_solution;
run;
data newt;
set newt;
rename COL1=x1 col2=x2;
run;
data params(type=est);
set newt;
_TYPE_='PARMS';
run;
proc nlp tech=NEWRAP gradcheck=detail gconv2=1e-21 inest=params;
min f;
parms x1, x2;
f = rosenbrock(x1,x2);
run;
Here's the whole code!
