Please help in below facility location code

[AAZ]

Hi Folks, 

I am new to Sas trying to solve below simple facility location, Please advice what I am missing,

Data	Wgt	X	Y
Boston	425	100	80
Providence	320	86	40
Springfield	220	20	60

Proc optmodel;

set city = {'BO','PR','SP'};
set Geo = {'X','Y'};

number Cor{city,Geo} =[100 80 86 40 20 60];

number Wgt{City}=[425 320 220];

var x{city,Geo} integer >= 0;
var y{city,Geo} integer >= 0;

Minimize MinWgtDis = sum{i in City}sum {j in Geo}
(sqrt((cor[i,j] - x[i,j])^2 + (cor[i,j] - y[i,j])^2 ))* Wgt[i];
/*
con x{i in city,j in Geo}: X[i,j]<=200;
con y{i in city,j in Geo}: y[i,j]<=200;*/

solve;

print MinWgtDis x y;
quit;

[RobPratt]

Here's one way to solve it as an unconstrained problem with nonlinear objective and two variables:

data indata;
   input city $12. Wgt X Y;
   datalines;
Boston      425 100  80
Providence  320  86  40
Springfield 220  20  60
;

proc sgplot data=indata;
   bubble x=x y=y size=wgt / datalabel=city;
run;

proc optmodel;
   set <str> CITIES;
   num wgt {CITIES};
   num xc {CITIES};
   num yc {CITIES};
   read data indata into CITIES=[city] wgt xc=x yc=y;

   var X, Y;

   min MinWgtDis = sum {i in CITIES} (sqrt((xc[i] - X)^2 + (yc[i] - Y)^2)) * wgt[i];

   solve;
   print MinWgtDis X Y;
   create data sganno from function='text' drawspace='datavalue' label='facility' x1=X y1=Y;
quit;

proc sgplot data=indata sganno=sganno;
   bubble x=x y=y size=wgt / datalabel=city;
run;

The resulting plot provides a sanity check on the optimal solution returned.

[AAZ]

Wow, that's super! Thanks

I am trying to solve it without defining data table, started to get some numbers but am still missing something, can you help me figuring out the error. I think coz I don't have a seperate variable for Xi and Yi it's not working right

Proc optmodel;

set city = {'BO','PR','SP'};
set Geo = {'Xi','Yi'};

number Cor{city,Geo} =[100 80 86 40 20 60];
number Wgt{city}=[425 320 220];

var x >= 0 ;
var y >= 0;

Minimize MinWgtDis = sum{i in city} sum {j in Geo}
((sqrt((cor[i,j] - x)^2 + (cor[i,j] - y)^2 ))* Wgt[i]);

solve;

print MinWgtDis x y ;

quit;

[RobPratt · Accepted solution]

You can use the EXPAND statement to see how the objective function you defined gets populated:

Var x >= 0                                                                                
Var y >= 0                                                                                
Minimize MinWgtDis=SQRT((-x + 100)**2 + (-y + 100)**2)*425 + SQRT((-x + 80)**2 + (-y + 80)
**2)*425 + SQRT((-x + 86)**2 + (-y + 86)**2)*320 + SQRT((-x + 40)**2 + (-y + 40)**2)*320  
 + SQRT((-x + 20)**2 + (-y + 20)**2)*220 + SQRT((-x + 60)**2 + (-y + 60)**2)*220          

You can see that the formula for MinWgtDis is not really what you intended.  Here's a modification of your code that instead does what you want:

Proc optmodel;
set city = {'BO','PR','SP'};
set Geo = {'Xi','Yi'};
number Cor{city,Geo} =[100 80 86 40 20 60];
number Wgt{city}=[425 320 220];
var x {Geo} >= 0 ;

Minimize MinWgtDis = sum{i in city} sqrt(sum {j in Geo} (cor[i,j] - x[j])^2)* Wgt[i];
 
solve;
print MinWgtDis x;
quit;

[AAZ]

how to use the expand statment to see how the function is populated ?

[RobPratt]

What I showed above came from the simplest invocation, which expands the whole model:

expand;

But you can also expand just the objective as follows:

expand MinWgtDis;

The full documentation for the EXPAND statement is here:

http://go.documentation.sas.com/?docsetId=ormpug&docsetVersion=14.2&docsetTarget=ormpug_optmodel_syn...

[AAZ]

I actually added two addional lines, but when I write expand, it expands only the optmization formula 

number SumCOG = sum(of Wgt[*]);

number xcog{geo} = sum{i in city} (sum{j in Geo} (Wgt[i]/SumCOG)*cor[i,j]);
expand ;

XCOG is not calculated right and want  trace it, but it doesn't show up in the expand result.

[RobPratt]

EXPAND does not apply here, but you can use the PRINT statement to inspect parts of your expression.  For example,

print SumCOG;

returns the correct value 965.

The issue with your NUM declaration is that you have used Geo both as an index set on the left-hand side and as a summation index set on the right-hand side.  I think what you intended is this:

num xcog{j in geo} = sum{i in city} ((Wgt[i]/SumCOG)*cor[i,j]);
print xcog;

SAS Output

[1]	xcog
Xi	77.119
Yi	62.176

By the way, the value of MinWgtDis is higher than the minimum value returned by the solver, which you can confirm as follows:

for {j in Geo} x[j] = xcog[j];
print MinWgtDis;

[AAZ]

that's extremly helpful Rob,I need to work on the summation notation. 

when you added 

"for {j in Geo} x[j] = xcog[j];
print MinWgtDis;"

 it means untill j holds a value, go back and use the new xcog instead of x  to calculate the minimum cost, right ?

and to print each and every indvidual X,Y, not the summed one.

I have to put this"sqrt(sum {j in Geo} (cor[i,j] - x[j])^2)* Wgt[i];"  in  avariable and then print it

[RobPratt]

The FOR loop sets x[j] to xcog[j] for each j in Geo.  The subsequent PRINT statement prints the new value of MinWgtDis, which gets automatically updated whenever any part of its definition gets updated.

A good source for learning about PROC OPTMODEL is SAS/OR 14.2 User's Guide: Mathematical Programming Examples.

[AAZ]

Thank you rob

[AAZ]

Hi Rob,

I added this code to your solution, it works fine, is there is a better way to do ?

/*Calculating CoG Coordinates */
num SumCOG = sum(of Wgt[*]);
for {i in CITIES} X = sum{i in CITIES} (Wgt[i]*xc[i])/SumCOG ;
for {i in CITIES} Y = sum{i in CITIES} (Wgt[i]*yc[i])/SumCOG ;

print MinWgtDis SumCOG X Y;
expand;

[RobPratt]

The FOR loop over CITIES is unnecessary, and you are performing the exact same calculation multiple times.  Instead, you can just do this:

X = sum{i in CITIES} (Wgt[i]*xc[i])/SumCOG ;
Y = sum{i in CITIES} (Wgt[i]*yc[i])/SumCOG ;