MILP Infeasibility

[Santha]

Hi 

I have this code below. trying to pick "p" number of facilities from a list of facilities. this works perfectly fine when I have p as 70 but when I change it to 60 it returns infeasibility. I want to identify which constraint/customer is causing it so that I can make necessary changes. I tried IIS=TRUE but it did not help. Earlier SAS used to give me which customer was causing infeasibility and i added fixes to it. Now I am not sure what causes infeasibility. Any help?

proc optmodel;
set DIMS=1..2;
set CUSTOMERS;
set FACILITIES;
num a {CUSTOMERS,DIMS};
num demand {CUSTOMERS}; 

num f {FACILITIES,DIMS};
num SiteCapacity {FACILITIES};
str FacilitySiteName {FACILITIES};

num p=70; /*Set how many Facilities you want to be OPENED*/

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


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

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

**print a;
**print SiteState;
**print FacilitySiteName;
**print demand;

/*str SiteRegion{CUSTOMERS};
read data STDOPT.COGMODELINGDATA intP CUSTOMERS=
[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> SiteRegion ;

str SiteComment1{CUSTOMERS};
read data STDOPT.COGMODELINGDATA into CUSTOMERS=
[_N_] {d in DIMS} <a[_N_,d]=col('a'||d)> SiteComment1;*/


/* 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],'M');

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

   min Cost
       = sum {<i,j> in CUSTOMERS_FACILITIES} dist[i,j] * Assign[i,j]*demand[i];
   
   /* 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 &SiteCapacity 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 the MILP */
 solve with milp /timetype=real;

This is the error message as well 

solve with milp /timetype=real;
NOTE: Problem generation will use 4 threads.
NOTE: The problem has 7341 variables (0 free, 0 fixed).
NOTE: The problem has 7341 binary and 0 integer variables.
NOTE: The problem has 7582 linear constraints (6937 LE, 645 EQ, 0 GE, 0 range).
NOTE: The problem has 21213 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The initial MILP heuristics are applied.
NOTE: The MILP presolver value AUTOMATIC is applied.
NOTE: The MILP presolver removed 90 variables and 91 constraints.
NOTE: The MILP presolver removed 199 constraint coefficients.
NOTE: The MILP presolver modified 0 constraint coefficients.
NOTE: The presolved problem has 7251 variables, 7491 constraints, and 21014 constraint coefficients.
NOTE: The MILP solver is called.
NOTE: The parallel Branch and Cut algorithm is used.
NOTE: The Branch and Cut algorithm is using up to 4 threads.
          Node   Active   Sols    BestInteger      BestBound      Gap    Time
             0        1      0              .         207564        .       0
             0        1      0              .         207564        .       0
NOTE: Infeasible.

[RobPratt]

Can you please share the data sets?

[Santha]

HI Rob

You think this could be in data problem?

[RobPratt]

I was asking for the data so that I can run your code and give you a specific diagnosis of infeasibility.

[Santha]

ok. let me share the data sets shortly. thank you

[Santha]

Rob

Here you go. This is a reduced data set on the cog customers. I ran the model on this and still got infeasibility

[RobPratt]

I get an optimal solution with that code and data:

NOTE: Problem generation will use 4 threads.
NOTE: The problem has 6251 variables (0 free, 0 fixed).
NOTE: The problem has 6251 binary and 0 integer variables.
NOTE: The problem has 6373 linear constraints (5845 LE, 528 EQ, 0 GE, 0 range).
NOTE: The problem has 17939 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The OPTMODEL presolver is disabled for linear problems.
NOTE: The initial MILP heuristics are applied.
NOTE: The MILP presolver value AUTOMATIC is applied.
NOTE: The MILP presolver removed 135 variables and 97 constraints.
NOTE: The MILP presolver removed 258 constraint coefficients.
NOTE: The MILP presolver modified 0 constraint coefficients.
NOTE: The presolved problem has 6116 variables, 6276 constraints, and 17681 constraint
      coefficients.
NOTE: The MILP solver is called.
NOTE: The parallel Branch and Cut algorithm is used.
NOTE: The Branch and Cut algorithm is using up to 4 threads.
          Node   Active   Sols    BestInteger      BestBound      Gap    Time
             0        1      0              .       18457359        .       0
             0        1      1       18570314       18457359    0.61%       0
             0        1      2       18510332       18510332    0.00%       0
             0        0      2       18510332       18510332    0.00%       0
NOTE: The MILP solver added 7 cuts with 346 cut coefficients at the root.
NOTE: Optimal.
NOTE: Objective = 18510331.726.
284  quit;
NOTE: PROCEDURE OPTMODEL used (Total process time):
      real time           1.10 seconds
      cpu time            1.15 seconds

[Santha]

can you try with p=60 ? and 

set CUSTOMERS_FACILITIES = {i in CUSTOMERS, j in FACILITIES: dist[i,j] <= 120};

this is what i got infeasible with this data set

[RobPratt]

Yes, I get infeasible with p = 60.  In fact, even the LP relaxation is infeasible, as you can see by calling SOLVE RELAXINT.  Both LP IIS and MILP IIS do not return a particularly small IIS in this case, but here is an alternative approach that you might find helpful.  Omit the Have_This_Many_FAC_OPEN constraint and instead minimize the number of facilities built:

*   con Have_This_Many_FAC_OPEN:
      Sum{j in FACILITIES} Build[j] <= p;

   min NumBuild = Sum{j in FACILITIES} Build[j];
   solve;

The MILP solver returns an optimal objective value of 61:

NOTE: Problem generation will use 4 threads.
NOTE: The problem has 6251 variables (0 free, 0 fixed).
NOTE: The problem has 6251 binary and 0 integer variables.
NOTE: The problem has 6372 linear constraints (5844 LE, 528 EQ, 0 GE, 0 range).
NOTE: The problem has 17532 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The OPTMODEL presolver is disabled for linear problems.
NOTE: The initial MILP heuristics are applied.
NOTE: The MILP presolver value AUTOMATIC is applied.
NOTE: The MILP presolver removed 181 variables and 142 constraints.
NOTE: The MILP presolver removed 304 constraint coefficients.
NOTE: The MILP presolver modified 0 constraint coefficients.
NOTE: The presolved problem has 6070 variables, 6230 constraints, and 17228 constraint
      coefficients.
NOTE: The MILP solver is called.
NOTE: The parallel Branch and Cut algorithm is used.
NOTE: The Branch and Cut algorithm is using up to 4 threads.
          Node   Active   Sols    BestInteger      BestBound      Gap    Time
             0        1      4     81.0000000     60.5000000   33.88%       0
NOTE: The MILP solver's symmetry detection found 4527 orbits. The largest orbit contains 28
      variables.
             0        1      8     61.0000000     60.5000000    0.83%       0
             0        0      8     61.0000000     61.0000000    0.00%       0
NOTE: Optimal.
NOTE: Objective = 61.

So any value of p < 61 will make the problem infeasible.

Another idea is to introduce explicit slack variables and minimize the sum of slacks.  Keep all of the original constraints, except replace the assign_def constraints as follows:

   var Unassigned {CUSTOMERS} >= 0;
   con assign_def {i in CUSTOMERS}:
      sum {<(i),j> in CUSTOMERS_FACILITIES} Assign[i,j] + Unassigned[i] >= 1;
   min NumUnassignedCustomers = sum {i in CUSTOMERS} Unassigned[i];

   solve;
   print {i in CUSTOMERS: Unassigned[i].sol > 0.5} Unassigned;

The resulting optimal solution has optimal objective value 1, with Unassigned[73] = 1, which suggests changing something about customer 73 if you want the problem to be feasible.

[Santha]

Rob

Ok . So we have two objectives? one is to minimize the build facilities. and another one is cost. Right?

[Santha]

I like the first approach of minimizing the facilities. This implies that the cost will also be reduced. Pls correct me if I am wrong

[RobPratt · Accepted solution]

I have introduced two new objectives beyond your original objective, but the solver uses only one at a time.  By default, it uses the one that was declared last, but you can also use the OBJ option in the SOLVE statement to explicitly specify which objective to use.

[Santha]

Thank you Rob very much. Appreciate your help as always

[Santha]

I am little confused. Now, we  have two objective functions.

min Cost= sum {<i,j> in CUSTOMERS_FACILITIES} dist[i,j] * Assign[i,j]*demand[i];
min NumBuild = Sum{j in FACILITIES} Build[j];

Since numbuild is the last , SAS solves for it. But if i Specify solve OBJECTIVE Cost, it would solve for the cost objective.

What I want the model to do is to minimize both Cost and NumBuild really. So, in that case, if I solve for the NumBuild Objective it will not be optimial cost one obviously. So, can i specify both  like this? 

solve with milp objective Cost,NumBuild /timetype=real;

[Santha]

Another option - Is it a good idea to run the min Numbuild objective function. Take that number and run the model for min cost? 

If yes, how do we do it?