Using Proc Optmodel for Set Covering Problem

[mhflah]

Hi all,

I am working on a set covering project where I am trying to locate facilities within a 25-minute distance from demand points. All destinations within 25-minutes of an origin were generated in GIS software external to SAS. 

So really, right now what I have is: Facility IDs, Destination ID, Destination Population, Capacity.

Simplified version would look like:

data have;

   input originid destinationid pop10 capacity;

   datalines;

1 1 15 22

1 2 5 22

2 2 5 22

2 3 17 22

2 4 4 22

3 3 17 22

3 4 4 22

4 2 5 22

4 4 4 22

;

What I want would be for origin 1 and 3 to be selected since that covers all points 1-4 without exceeding the capacity of 22, while excluding origin id 2, since placement there would exceed 22, and 4, since 1 and 3 serve more while still covering each demand point. Overlap is okay. 

At the end of the day, I'm starting to think this may be a data prep problem on my end. This is how I've been conceptualizing it, but I'm not sure if my have data set is telling the same story to the computer:

Thank you for your time. 

[RobPratt · Accepted solution]

Please see if the following does what you want:

proc optmodel;
   set <num,num> ORIGINS_DESTINATIONS;
   read data have into ORIGINS_DESTINATIONS=[originid destinationid];

   set ORIGINS;
   num capacity {ORIGINS};
   read data have into ORIGINS=[originid] capacity;

   set DESTINATIONS;
   num demand {DESTINATIONS};
   read data have into DESTINATIONS=[destinationid] demand=pop10;

   var IsOpen {ORIGINS} binary;
   var Serve {ORIGINS_DESTINATIONS} binary;

   con Cover {d in DESTINATIONS}:
      sum {<o,(d)> in ORIGINS_DESTINATIONS} Serve[o,d] >= 1;
   con CapacityCon {o in ORIGINS}:
      sum {<(o),d> in ORIGINS_DESTINATIONS} demand[d] * Serve[o,d] <= capacity[o] * IsOpen[o];

   /* primary objective: minimize number of open origins */
   min NumOpen = sum {o in ORIGINS} IsOpen[o];

   /* secondary objective: maximize total coverage of demands */
   max TotalCoverage = sum {<o,d> in ORIGINS_DESTINATIONS} demand[d] * Serve[o,d];

   solve obj NumOpen;
   print NumOpen TotalCoverage;
   print Serve;
   print IsOpen {o in ORIGINS} (sum {<(o),d> in ORIGINS_DESTINATIONS} demand[d] * Serve[o,d]);

   num minNumOpen;
   minNumOpen = NumOpen.sol;
   con ObjectiveCut:
      NumOpen <= minNumOpen;

   solve obj TotalCoverage;
   print NumOpen TotalCoverage;
   print Serve;
   print IsOpen {o in ORIGINS} (sum {<(o),d> in ORIGINS_DESTINATIONS} demand[d] * Serve[o,d]);
quit;

NumOpen	TotalCoverage
2	41

Serve
	1	2	3	4
1	1	1
2		0	1	1
3			0	0
4		0		0

[1]	IsOpen
1	1	20
2	1	21
3	0	0
4	0	0

[mhflah]

This is exactly what I was trying to do.

Thank you, Rob!