Here is code to solve the multisource Weber problem, where m is the number of facilities. It is based on the code I posted in the other thread, but you can modify as needed for your data and distance function. Because the multisource problem is nonconvex, I used the multistart option for the NLP solver and also added optional code to post-process the solution to make further local improvements.
%let m = 3;
proc optmodel;
set DIMS = 1..2;
set CUSTOMERS;
set FACILITIES = 1..&m;
num a {CUSTOMERS, DIMS};
num demand {CUSTOMERS};
read data cdata into CUSTOMERS=[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> demand;
num Xlb {d in DIMS} = min {i in CUSTOMERS} a[i,d];
num Xub {d in DIMS} = max {i in CUSTOMERS} a[i,d];
var X {FACILITIES, d in DIMS} >= Xlb[d] <= Xub[d];
var W {CUSTOMERS, FACILITIES} >= 0;
impvar Distance {i in CUSTOMERS, j in FACILITIES} =
sqrt(sum {d in DIMS}(a[i,d]-X[j,d])^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;
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;
create data Xdata from [j]=FACILITIES {d in DIMS} <col('X'||d)=X[j,d].sol>;
create data assignments from [j i]={j in FACILITIES, i in CUSTOMERS_j[j]} W[i,j]
x1=a[i,1] y1=a[i,2] x2=X[j,1] y2=X[j,2];
quit;
proc sgplot data=assignments;
vector x=x2 y=y2 / xorigin=x1 yorigin=y1 group=j;
run;
... View more
