Looking for Literature for Center of Gravity problem

[Santha]

Hi.

I am trying to see if I can solve the center of gravity location problem using SAS. So I am looking for articles / literature to set my model. 

As an overview, I have a bunch of nodes , each of them with its own demand (weight). I have latitudes and longitudes of them. 

I am looking to find a center of gravity point that will minimize the "weighted-distance" from this center of gravity point to all nodes. 

And as an extension, want to see if i can choose from a list of candidates that will minimize the weighted-distance.

Any help is appreciated

[RobPratt · Accepted solution]

This is known as the Weber problem and can be solved as an unconstrained convex nonlinear optimization problem.  See this SAS Global Forum 2008 paper, which also discusses the (nonconvex) multisource Weber problem.

[Santha]

Thanks Rob. I will look into it. 

[Santha]

Rob 

I did see that. It gave me a good insight. Thank you. 

I was going to ask you for any example models in SAS related to this that I can learn so that I could better formulate it.

Thanks again

[RobPratt]

I happened to already have code for this:

proc optmodel;
   set DIMS = 1..2;
   set CUSTOMERS;
   num a {CUSTOMERS, DIMS};
   num demand {CUSTOMERS};
   read data indata into CUSTOMERS=[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> demand;

   var X {DIMS};
   min Z = sum {i in CUSTOMERS} demand[i]*sqrt(sum {d in DIMS}(a[i,d]-X[d])^2);

   solve;
   print X;

   create data sganno from [i]=CUSTOMERS x1=a[i,1] y1=a[i,2] x2=X[1] y2=X[2]
      drawspace='datavalue' function='line';
quit;

proc sgplot data=cdata sganno=sganno;
   scatter x=a1 y=a2;
run;

[Santha]

Rob that is really amazing.thank you

[Santha]

Rob. 

I went into your model and also did some research to understand basic concepts about these types of problems. 

Now I am ready to deploy SAS and test it on my small model. I have two columns called Sitename (that are nothing but customer names ) and Demand. I have 6 Sitenames and demand for each 6 of them. My objectives is to pick "n" points that will minimize the sum of the weighted distance to all the customers. Once this is set up correctly, I shall go on with choosing from a list of candidate facilities. 

Sitename	SiteParameter1
DEN	800
SFO	150
EWR	300
SEA	200
MIA	50
MCI	500
ORD	1000

I will have a set of customers (set CUSTOMERS;)

I will have a numeric parameter called demand that is indexed over CUSTOMERS  (num demand {CUSTOMERS};)

The other numeric parameter that I see in your code num a{CUSTOMERS,DIMS} and your read data I am not quite sure I follow it. I am understanding that to be another parameter which will be like Customer1, Customer2 etc? 

read data indata into CUSTOMERS=[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> demand;  

Can you help me how to read data / understand read data ? Also, model will give the optimum X,so that will be one of the 7 sites.  Correct?

[RobPratt]

In the example code I provided, the indata data set has two columns (a1 and a2) that correspond to the location coordinates of each customer.  For you, these would be latitude and longitude, and you would probably want to use the GEODIST function to compute distance instead of the SQRT formula that I used to compute Euclidean distance.  Solving the optimization model yields the coordinates (X[1] and X[2]) for a new location that minimizes the demand-weighted sum of distances from customers to this new location.  It will not be the same location as an existing customer unless that customer happens to be located at the center of gravity.

[Santha]

Thank you Rob for the insights. I am clear now. I will set up the model and let you know. 

This is very interesting. Thank you again.

[Santha]

Rob

I wanted to let you know that I was able to adapt the code to my needs. and was able to get what I exactly wanted. 

Thank you very much.

[RobPratt]

Glad to hear that, Santha!

[Santha]

I have a full model with 11,000 customers. Its been like 45 minutes and is still running.

I want model to pick 1 answer. 

How do I set the tolerance / relative gap so that I can live with that % . Also, a time limit until the model can run like an hour or so and give the answers after 1 hr mark.?

[Santha]

Trying FEASTOL and MAXTIME based on NLP options. Any other ideas?

[RobPratt]

Yes, MAXTIME=3600 will stop the solver after one hour.  You can also use FEASTOL and/or OPTTOL.  Are you able to share code and data?

[Santha]

Rob - here is the code. 

%let p=2; 
proc optmodel;

set DIMS=1..2;
set CUSTOMERS;
set FACILITIES = 1..&p;
num a {CUSTOMERS,DIMS};
num demand {CUSTOMERS};

read data STDOPT.COGModelingData into CUSTOMERS=
[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> demand;

print demand;
print a;

num Xlb {d in DIMS} = min {i in CUSTOMERS} a[i,d];
num Xub {d in DIMS} = max {i in CUSTOMERS} a[i,d];
**print Xlb; **print Xub;

var X {FACILITIES, d in DIMS} >= Xlb[d] <= Xub[d];
var W {CUSTOMERS, FACILITIES} >= 0;

impvar Distance {i in CUSTOMERS, j in FACILITIES} = sum{d in DIMS} GEODIST(a[i,1],a[i,2],X[j,1],X[j,2]);

min Z = sum {i in CUSTOMERS, j in FACILITIES} W[i,j]*Distance[i,j];
con DemandCon {i in CUSTOMERS}: sum {j in FACILITIES} W[i,j] = demand[i];

solve with nlp / ms FEASTOL=1e-1;
print X;
put _OBJ_=;

/* post-processing: assign each customer to closest facility */
set CUSTOMERS_j {FACILITIES} init {};
num minDistance, argminDistance;
for {i in CUSTOMERS} do;
minDistance = constant('BIG');
argminDistance = .;
for {j in FACILITIES} do;
if minDistance > Distance[i,j] then do;
minDistance = Distance[i,j];
argminDistance = j;
end;
end;
for {j in FACILITIES} W[i,j] = 0;
W[i,argminDistance] = demand[i];
CUSTOMERS_j[argminDistance] = CUSTOMERS_j[argminDistance] union {i};
end;
put _OBJ_=;

/* post-processing: solve each facility separately */
min SingleFacilityObjective {j in FACILITIES} = sum {i in CUSTOMERS_j[j]} demand[i]*Distance[i,j];
problem SingleFacilityProblem {j in FACILITIES} include
{d in DIMS} X[j,d]
SingleFacilityObjective[j];
for {j in FACILITIES} do;
put j=;
use problem SingleFacilityProblem[j];
solve;
end;
put (sum {j in FACILITIES} SingleFacilityObjective[j])=;
print X;