Here are sample data and code for the p-median problem. If you want to force certain facilities to be open, you can use the FIX statement to fix the corresponding Buiild[j] variables to 1. Note that this formulation is linear because the nonlinear distances are absorbed into the objective coefficients.
%let p = 5;
%let NumCustomers = 200;
%let xmax = 100;
%let ymax = 100;
/* generate random customer locations */
data cdata;
call streaminit(1);
do i = 1 to &NumCustomers;
x = rand('UNIF') * &xmax;
y = rand('UNIF') * &ymax;
weight = 1;
output;
end;
run;
/* plot customer locations */
proc sgplot data=cdata noautolegend aspect=1;
xaxis display=(nolabel);
yaxis display=(nolabel);
scatter y=y x=x;
run;
/* solve the problem by using MILP */
proc optmodel;
set CUSTOMERS;
/* x and y coordinates of customers */
num x {CUSTOMERS};
num y {CUSTOMERS};
num weight {CUSTOMERS};
read data cdata into CUSTOMERS=[i] x y weight;
/* distance between customers i and j */
num distance {i in CUSTOMERS, j in CUSTOMERS}
= sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2);
/* declare decision variables */
var Build {CUSTOMERS} binary;
var Assign {CUSTOMERS, CUSTOMERS} binary;
/* declare objective */
min TotalDistance
= sum {i in CUSTOMERS, j in CUSTOMERS} weight[i] * distance[i,j] * Assign[i,j];
/* build exactly p facilities */
con Cardinality:
sum {j in CUSTOMERS} Build[j] = &p;
/* each customer assigned to exactly one site */
con AssignOnce {i in CUSTOMERS}:
sum {j in CUSTOMERS} Assign[i,j] = 1;
/* if customer i assigned to site j, then facility must be built at j */
con AssignImpliesBuild {i in CUSTOMERS, j in CUSTOMERS}:
Assign[i,j] <= Build[j];
* con AssignImpliesBuild {j in CUSTOMERS}:
sum {i in CUSTOMERS} Assign[i,j] <= card(CUSTOMERS) * Build[j];
/* solve the MILP */
solve;
/* create a data set for use by PROC SGPLOT */
create data solution from
[customer site]={i in CUSTOMERS, j in CUSTOMERS: Assign[i,j].sol > 0.5}
x1=x[i] y1=y[i] x2=x[j] y2=y[j];
quit;
/* plot optimal solution */
proc sgplot data=solution noautolegend aspect=1;
xaxis display=(nolabel);
yaxis display=(nolabel);
scatter x=x1 y=y1;
vector x=x2 y=y2 / xorigin=x1 yorigin=y1 noarrowheads;
run;