p-Median Weber

[Santha]

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?

[RobPratt · Accepted solution]

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;