topic Re: MILP Infeasibility in Mathematical Optimization, Discrete-Event Simulation, and OR

MILP Infeasibility

Santha — Thu, 05 Aug 2021 20:20:42 GMT

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.

Re: MILP Infeasibility

RobPratt — Thu, 05 Aug 2021 20:23:22 GMT

Can you please share the data sets?

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 20:30:31 GMT

HI Rob

You think this could be in data problem?

Re: MILP Infeasibility

RobPratt — Thu, 05 Aug 2021 20:37:09 GMT

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

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 20:37:47 GMT

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

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 20:51:30 GMT

Rob

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

Re: MILP Infeasibility

RobPratt — Thu, 05 Aug 2021 21:10:22 GMT

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

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 21:18:10 GMT

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

Re: MILP Infeasibility

RobPratt — Thu, 05 Aug 2021 21:49:01 GMT

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.

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 21:43:39 GMT

Rob

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

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 21:52:00 GMT

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

Re: MILP Infeasibility

RobPratt — Thu, 05 Aug 2021 21:52:02 GMT

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.

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 23:10:53 GMT

Thank you Rob very much. Appreciate your help as always

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 23:46:11 GMT

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;

Re: MILP Infeasibility

Santha — Thu, 05 Aug 2021 23:48:15 GMT

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?

Re: MILP Infeasibility

RobPratt — Fri, 06 Aug 2021 00:44:36 GMT

No, the MILP solver optimizes only one objective at a time. The black-box solver can handle multiple objectives, but with no guarantee of optimality. See this doc example: SAS Help Center: Multiobjective Optimization

Re: MILP Infeasibility

RobPratt — Fri, 06 Aug 2021 00:46:58 GMT

Yes, that is a natural approach when one objective is more important than the other. Here's how you can do it:

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=60; /*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 ;

   /* 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;

   /* 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];

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

   num p;
   p = NumBuild;
   con Have_This_Many_FAC_OPEN:
      NumBuild <= p;

   min Cost
      = sum {<i,j> in CUSTOMERS_FACILITIES} dist[i,j] * Assign[i,j]*demand[i];
   solve with milp / primalin;
quit;

Re: MILP Infeasibility

Santha — Fri, 06 Aug 2021 01:25:28 GMT

Rob

this is awesome. I was so glad that i was also on the same path coding. It worked perfectly fine the way I wanted to.

Thank you. You are simply amazing

Re: MILP Infeasibility

Santha — Fri, 27 Aug 2021 20:34:26 GMT

Hi Rob

For the same model in this thread, I tried to do this and got "Out of memory" error

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

what i want to do is for the model to pick exactly 2 facilities from a pool of candidate facilities.

So in order to do that I set this radius as like a big number like 3000 miles so that there model would solve. It gave me out of memory error. I understand that as giving a radius of 3000 miles means almost all possible combinations of faciliy- customer points. I tried to reduce that number to a point where it is like 550 or 600 then the model would not give me this "out of memory" error. But it would say infeasible which is understandable. So, how do i set this radius which I am using to pre_pick this CUSOMER_FACILITIES list within a certain range and yet solve this p number of facilities problem

Re: MILP Infeasibility

RobPratt — Fri, 27 Aug 2021 21:10:36 GMT

You can use the following statement to determine the smallest feasible distance threshold:

   num threshold = min {j1 in FACILITIES, j2 in FACILITIES: j1 < j2} max {i in CUSTOMERS} min(dist[i,j1], dist[i,j2]);

This threshold value turns out to be 1018.4387857. Now use threshold in the SET statement:

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

The resulting solution has optimal objective value 375267751.16.