defining inequality constraint for PSO

[CerditoSalvaje]

aloa everybody,

I'm trying to code PSO in SAS and I have problems with the code defining constraints inside the objective function (using the penalty method). My algorithm works well on unconstrained problem, but when I add a constraint it finds a wedge solution...

The theoritical solution to my problem are : (x1*,x2*)=(3.75;5.625).

proc iml;

/***parameters and variables***/
n=20;						/*number of particules*/
p=2;						/*number of dimensions*/
wmax=0.9;					/*inertia weight*/
wmin=0.4;					/*inertia weigth*/
c1=2.05;						/*acceleration factor*/
c2=2.05;						/*acceleration factor*/
phi=c1+c2;
chi= 2/abs(2-phi-sqrt(phi**2 -4*phi));				/*constriction factor*/
maxt=1000;					/*maximum iteration*/
maxrun=2;					/*trials of the algorithm*/
X=j(n,p+1,0);				/*initialization of the particle's matrix*/


start f(x);
con=100*x[1]+200*x[2]-1500;
y= 3*(x[1]**0.25)*(x[2]**0.75);
pen=10**9;
penalty=0;
if con>0 then penalty=pen*con;
else e=0;
fc=y+penalty;
return(fc);
finish f;


xmin=0;
xmax=100;
lb=j(1,p,xmin);
ub=j(1,p,xmax);

/*other useful matrixes*/
ffmin=j(maxt,maxrun,0);
fft=j(1,maxrun,0);
fff=j(1, maxrun,0);
rgbest=j(maxrun,p,0);
nb_t=j(2, maxrun,0);

do run=1 to maxrun;
	nb_t[1,run]=run;
					/***initialization of the pso***/
	do i=1 to n;
		do j=1 to p;
	/*sample in a uniform distribution to initialize the particles position*/
			call randgen(u,'uniform');
			X[i,j]=lb[j]+u*(ub[j]-lb[j]);
		end;
	/*to fill the third column of X (which corresponds to the output of the function), we evaluate the fitting for each particule*/
		X[i,p+1]=f(X[i,]);
	end;

	pbest=X;			/*initialization of pbest so we can iterate then*/
	V=0.1*X;			/*initialization of the velocity (10%) */

	/*we find the current minimum of f, given the current values of X (if we want to minimize the function)*/
	titi=min(pbest[,p+1]);
	/*then we iterate to assign the coordinates of this minimum to gbest*/
	do i=1 to n;
		if pbest[i,p+1]=titi then gbest=pbest[i,];
		else e=0;
	end;
	*print X, pbest, gbest;
/**********************************/
	t=1;							/*resetting the iteration counter*/
	tolerance=1;					/*resetting the tolerance(stopping criteria)*/
	do until(tolerance < (10**(-12))|t=maxt);

		w=wmax-(wmax-wmin)*t/maxt;			/*updating the inertial weight parameter*/

		/*for all the particles in all dimensions : */
		do i=1 to n;
			do j=1 to p;

				/*velocity updates*/
				call randgen(u1,'uniform');
				call randgen(u2,'uniform');
				V[i,j]=chi*(w*V[i,j]+c1*u1*(pbest[i,j]-X[i,j])+c2*u2*(gbest[1,j]-X[i,j]));

				/*position updates*/
				X[i,j]= X[i,j]+V[i,j];

				/*handling boundary violations*/
				if X[i,j]<lb[j] then X[i,j]=lb[j];
				else e=0;
				if X[i,j]>ub[j] then X[i,j]=ub[j];
	 			else e=0;

			end;
		end;
		
		/*computing fitness for all the particules + updating pbest*/ 
		do i=1 to n;

			X[i,p+1]=f(X[i,]);

			if X[i,p+1]<pbest[i,p+1] then pbest[i,]=X[i,];
			else e=0;

		end;
		
		ffmin[t,run]=gbest[,p+1];
		fft[run]=t;

		/* updating gbest*/
		titi=min(pbest[,p+1]);
		do i=1 to n;
			if pbest[i,p+1]=titi then gbest=pbest[i,];
			else e=0;
		end;

		/*calculating tolerance*/
		if t>100 then tolerance=abs(ffmin[(t-100),run]-gbest[,p+1]);
		else e=0;
		
		t=t+1;

	end;
	*print gbest;
	nb_t[2,run]=t;					/*storing the number of iterations it took the algorithm to converge*/
	fff[run]=gbest[,p+1];
	rgbest[run,]=gbest[,1:p];
end;

/*displaying our results*/
fff2=fff`;
bestfun=mean(fff2);			/*taking the best run among all of the runs of the simulation*/
bx1=mean(rgbest[1]);
bx2=mean(rgbest[2]);

*print bestrun, bestfun, best_variables;

titi=nb_t`;
average=mean(titi[,2]);

best_variables= bx1||bx2;

toto=bestfun || bx1 || bx2 || average;			/*we concatenate the main results*/

best_results= TableCreate({"Bestfun" "bx1" "bx2" "average"}, toto);
call TablePrint(best_results) template="PSO_Results";

quit;

Best regards,

[Rick_SAS]

I don't know what a PSO is, nor a "wedge solution." Can you clarify what you are trying to do?  What is the purpose of all the code that has loops and random values?

I see a function 

y= 3*(x[1]**0.25)*(x[2]**0.75)

and your discription makes it sound like you want to subject it to the constraints

100*x[1]+200*x[2]-1500 <= 0?

I assume that you intend to also constrain
x[1] > 0
x[2] > 0

This is a nonlinear optimization problem with boundary and linear constraints, which you can solve by using one of the NLP functions in IML. For example, if you want to maximize the function subject to the constraints, you can run the following:

proc iml;
start func(x);
   y= 3*(x[1]**0.25)*(x[2]**0.75);
   return(y);
finish;

con = {
1e-6 1e-6 . . ,          /* x[1] > 0 and x[2] > 0 */
.    .    . . ,
100  200 -1 1500    /* 100*x[1]+200*x[2]-1500 < 0 */
};
x0 = {2 3};    /* initial guess */
opt = {1 0};   /* maximization, noprint */
call nlpnra(rc, xOpt, "func", x0) BLC=con opt=opt;
print rc, xOpt;

The program find the solution

xOpt
3.7499782	5.6250109

[sbxkoenk]

@Rick_SAS : I think they are trying to do particle swarm optimization.

There are some SAS Global Forum (proceedings) papers on that topic.

Koen

[Rick_SAS]

Before you proceed any further, please consider vectorizing your code. You can allocate a vector of random numbers with one call, then operate on the vectors. When you avoid looping over individual elements, your code goes much faster. It is also easier to develop, debug, and explain to others. If you don't know how to vectorize PROC IML code, here is a short tutorial.