mi

[Santha]

Rob

Here is my code that I followed from SPARSE modeling for MILP. I have 9 facilities to choose from. I have 1,056 customers. I want the CUSTOMER_FACILITIES to be the place where I can set my radius. Also, I want to be able to limit the number of facilities to "p" that i set initially.  I have a constraint (con Have_This_Many_FAC_OPEN) for that purpose . Is that correct? Right now I have this commented. When I ran my model with distance <=10,000 the model ran. However, i do see that each customer is assigned to all 9 different facilities, which means that the constraint "assign_def" is not set up correctly and/OR the CUSTOMER_FACILITIES is not set up right. If I reduce my radius to 800 for example because technically any customer can be served from a facility within 800 miles within continental US. But it results in in feasibility. Not sure what I am doing wrong. The SiteCapacity for all facility is set to be a very big number so that right now there is no in feasibility because of this. I can change sitecapacity later.

%let p=4; /*Set how many Facilities you want to be OPENED*/
proc optmodel;
set DIMS=1..2;
set CUSTOMERS;
set FACILITIES;
num a {CUSTOMERS,DIMS};
num demand {CUSTOMERS}; 
num f {FACILITIES,DIMS};
num SiteCapacity {FACILITIES};

read data STDOPT.COGcustomers into CUSTOMERS=
[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> demand;

read data STDOPT.cogExistingFacilities into FACILITIES=
[_N_] {d in DIMS} <f[_N_,d]=col('f'||d)> SiteCapacity;

/* distance from customer i to facility j */
num dist {i in CUSTOMERS, j in FACILITIES}
= GEODIST(a[i,1],a[i,2],f[j,1],f[j,2]);

set CUSTOMERS_FACILITIES = {i in CUSTOMERS, j in FACILITIES: dist[i,j] <= 10000};
var Assign {CUSTOMERS_FACILITIES} binary;
var Build {FACILITIES} binary;

min Cost
= sum {<i,j> in CUSTOMERS_FACILITIES} dist[i,j] * Assign[i,j];
   
/* each customer assigned to exactly one site */
 con assign_def {i in CUSTOMERS}:
 sum {<(i),j> in CUSTOMERS_FACILITIES} Assign[i,j] = 1;

 /* if customer i assigned to site j, then facility must be built at j */
   con link {<i,j> in CUSTOMERS_FACILITIES}:
      Assign[i,j] <= Build[j];
/* each site can handle at most Site Capacity demand */
con capacity {j in FACILITIES}:
sum {<i,(j)> in CUSTOMERS_FACILITIES} demand[i] * Assign[i,j] <= SiteCapacity[j] * Build[j];
 
**con  Have_This_Many_FAC_OPEN:  Sum{j in FACILITIES} Build[j]=p;

solve with milp/timetype=real;

[RobPratt]

Yes, the commented constraint is correct, except that you need &p instead of p unless you also declare p as a num:

num p = &p;

I don't immediately see anything wrong with the rest of the code, so maybe something is wrong with the data.  Please provide the input data sets so that I can run this myself to investigate further.

[Santha]

ok thanks rob. updated the variable p. here is my data set in csv format.

two files (COGCUSTOMERS) and COGEXISTINGFACILITIES.

[RobPratt · Accepted solution]

Note that GEODIST defaults to kilometers.  If you want miles, include 'M' as a fifth argument.

When I run with distance threshold 10000, I get one facility per customer, as expected from the assign_def constraints.

When I run with distance threshold 800, I get infeasible, as you did.  These notes in the log indicate several customers that are not within 800 miles of any facility:

NOTE: The constraint 'assign_def[478]' is empty and infeasible.
NOTE: The constraint 'assign_def[486]' is empty and infeasible.
NOTE: The constraint 'assign_def[488]' is empty and infeasible.
NOTE: The constraint 'assign_def[489]' is empty and infeasible.
NOTE: The constraint 'assign_def[490]' is empty and infeasible.
NOTE: The constraint 'assign_def[491]' is empty and infeasible.
NOTE: The constraint 'assign_def[492]' is empty and infeasible.
NOTE: The constraint 'assign_def[902]' is empty and infeasible.

 The minimum feasible threshold turns out to be 926.21776127 miles, the distance from customer 492 to facility 6.

[Santha]

Rob. 

Ok. I missed the "m" as fifth argument. ran the model for 1000 miles . But if u look at the attached output ( filter for any one customer, i see 9 possibile (all 9) facilities as options. Should i print for which one is the winner by displaying the winner facility using Assign variable? 

[Santha]

When i printed the assign variable I do see that one customer is assigned to only one facility. so that is what you mentioned.

Let me see if p number of facilites are open as well

[Santha]

I tested the number of facilities and it worked as well. All good. Just two questions on outputs

1) i want to see only those customer-facility pairs where the assign variable is 1.  I tried this below by throwing a filter of Build = 1 but it did not work.

create data STDOPT.COG_Output(where=(Build=1)) from [i j]={i in CUSTOMERS,j in FACILITIES} 
      Cust_Lat=a[i,1] Cust_Lon=a[i,2] Fac_Lat=f[j,1] Fac_Lon=f[j,2] Distance=dist[i,j]
Assign SiteState=SiteState[i];

 2) How do I display the Sitename of the facility (the actual name) instead of the Sitenumber. like "Charleston" as opposed to number 9. 

[Santha]

oops. sorry. I got it. I got the wrong variable in my create data step . When i changed to Assign, it worked 

[RobPratt]

Comparison of Assign = 1 is dangerous because of numerical tolerances.  It is safer to do Assign > 0.5 like this:

create data STDOPT.COG_Output(where=(Assign>0.5)) from [i j]={i in CUSTOMERS,j in FACILITIES} 

or this:

create data STDOPT.COG_Output from [i j]={i in CUSTOMERS,j in FACILITIES: Assign[i,j].sol > 0.5} 

If you want to include SiteName in the output, first declare and read it:

str SiteName {FACILITIES};
read data STDOPT.cogExistingFacilities into FACILITIES=
   [_N_] {d in DIMS} <f[_N_,d]=col('f'||d)> SiteCapacity SiteName;

Then add it to the CREATE DATA statement:

create data STDOPT.COG_Output from [i j]={i in CUSTOMERS,j in FACILITIES: Assign[i,j].sol > 0.5} 
   Cust_Lat=a[i,1] Cust_Lon=a[i,2] Fac_Lat=f[j,1] Fac_Lon=f[j,2] Distance=dist[i,j]
   Assign SiteName=SiteName[j];

Note the [j] index because SiteName depends on j only and not on both i and j.

[Santha]

Rob

Thank you. that makes total sense. All questions are answered perfectly. Thank you again