## p-Median Weber

Rob

I am trying to extend my model to do this p-median.

What I have so far is model picks "x" number of facilities in the map.

I have a very finite set of pre-determined facilities (no more than 15). From this, I want to pick "x" . and (b) from this 15, I have already picked lets say 2 facilities that I want the model to use. I want the model to give me just one more, there by giving a total number of facilities opened to be 3.

I went through the paper (https://support.sas.com/resources/papers/proceedings/pdfs/sgf2008/203-2008.pdf)

I saw that one of the way to do is to have the facilities as customer. I tried that but did not quite achieve that with my code. I also tried one more approach where I store the list of facilities in an array with their lat and lon and was trying to have a binary variable for the site to be open or close. But realized that the first part is NLP and so i cannot have integer variables. I am a little lost here. I am thinking to solve NLP first but I want to restrict the solve to only the list of candidate facilities and not the entire map. Then do post processing to fine tune the assignment.  But i am not sure. Can you help me with your thoughts?

1 ACCEPTED SOLUTION

Accepted Solutions

## Re: p-Median Weber

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;
``````

## Re: p-Median Weber

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;
``````

Discussion stats