7
$\begingroup$

Algebraic Formulation - Note: All parameter values are subject to change, I have just used the numbers you see here as place holders for the time being. Additionally, there are a few other equations but I am highly confident that they are not causing the infeasibility problems that I am having. \begin{align}\min\quad&\sum_{i=1}^I\sum_{j=1}^J\sum_{t=1}^T(2+5 L_{t}^{ij}) I_{t,post}^{ij}\\\text{s.t.}\quad&I_{t,post}^{ij} = \frac{I_{t-1,treated}^2} {I_{t-1,treated}^2 + a^2} S_{t-1}^{ij}\\\quad&I_{t,treated} = (1 - 0.5L_{t}^{ij})I_{t,post}^{ij}\\\quad&S_{t}^{ij} = S_{t}^{ij} + sj_{t-1,k=2}^{ij} - I_{t-1,treated}\\\quad&\sum_{t=1}^T\sum_{i=1}^I\sum_{j=1}^J2L_{t}^{ij} \le B \\\quad&0 \le I_{t,post}^{ij},I_{t,treated}^{ij},S_{t}^{ij} \\ \quad&0 \le I_{0,treated}^{ij} \,\ S_{0}^{ij} \,\ j_{0,k=2}^{ij}\\\quad&L_{t}^{ij}=\,\text{ binary decision variable}\end{align}

I am running a MINLP model of invasive species control and I would ideally want all of my variables (except the binary variable L) to be in the NonNegativeReals domain, however, when I set this specification on any of the variables it forces some of my constraints to become infeasible. I can run the model with all variables in the Real domain and everything works out okay but it provides me with negative numbers which is not realistic in regards to the problems physical interpretation.

Consider the case where all variables are set to real and we look at the first few solutions for Constraint 4: Susceptible Recruitment under Output 1, we can see that at (1,1,3) that the variable begins taking negative values, and the culprit in the constraint is model.inf_treated[i,j,t-1].

# Constraint 4: Susceptible recruitment
def susceptible_advance_rule(model, i, j, t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.susceptible[i, j, t ] == model.susceptible[i, j, t-1] + model.juvsurv * model.juvenille[i, j, t-1, Kmax] - model.inf_treated[i,j,t-1]
model.susceptible_advancement = Constraint(model.Iset, model.Jset, model.Tset, rule=susceptible_advance_rule)

Output 1: model.susceptible in the 'Reals' domain

susceptible : Size=20, Index=susceptible_index
        Key       : Lower : Value               : Upper : Fixed : Stale : Domain
        (1, 1, 0) :  None :                50.0 :  None :  True :  True :  Reals
        (1, 1, 1) :  None :                53.2 :  None : False : False :  Reals
        (1, 1, 2) :  None :   19.28555555556563 :  None : False : False :  Reals
        (1, 1, 3) :  None :  -32.15963392758387 :  None : False : False :  Reals
        (1, 1, 4) :  None :  -50.31374332256311 :  None : False : False :  Reals
        (1, 2, 0) :  None :                30.0 :  None :  True :  True :  Reals
        (1, 2, 1) :  None :  30.700000000000003 :  None : False : False :  Reals
        (1, 2, 2) :  None :  14.750000000000497 :  None : False : False :  Reals
        (1, 2, 3) :  None : -15.149263681592133 :  None : False : False :  Reals
        (1, 2, 4) :  None : -29.438417445000606 :  None : False : False :  Reals

In the case, I would want to try and restrict the domain of model.susceptible to the NonNegativeReals and upon doing so I would expect that all of my model variables become greater than 0, but I get infeasibility in my constraints as a result, and my binary variables become non-binary as a symptom. The same infeasibility comes about if I then turn every one of my variables to exist in the domain of the NonNegativeReals, but I think the infeasibility comes from the conflict that arises from Constraint 4 when I restrict model.susceptible to be positive, but then I subtract away model.inf_treated.

Output 2: model.susceptible in the 'NonNegativeReals' domain

 susceptible : Size=20, Index=susceptible_index
        Key       : Lower : Value                  : Upper : Fixed : Stale : Domain
        (1, 1, 0) :     0 :                   50.0 :  None :  True :  True : NonNegativeReals
        (1, 1, 1) :     0 :                   53.2 :  None : False : False : NonNegativeReals
        (1, 1, 2) :     0 :     63.222379594592645 :  None : False : False : NonNegativeReals
        (1, 1, 3) :     0 :      58.85143663732718 :  None : False : False : NonNegativeReals
        (1, 1, 4) :     0 : 1.4845965610247883e-08 :  None : False : False : NonNegativeReals
        (1, 2, 0) :     0 :                   30.0 :  None :  True :  True : NonNegativeReals
        (1, 2, 1) :     0 :     30.700000000000003 :  None : False : False : NonNegativeReals
        (1, 2, 2) :     0 :     14.750000001550784 :  None : False : False : NonNegativeReals
        (1, 2, 3) :     0 :      6.753462538302196 :  None : False : False : NonNegativeReals
        (1, 2, 4) :     0 :  3.620342424247316e-10 :  None : False : False : NonNegativeReals

Infeasible Constraints When Variable Boundaries Appied

#Constraint 9: Treated Infestation
def treatment_rule(model, i, j, t):
    if t == 0:
        return Constraint.Skip
    else:
        return model.inf_treated[i, j, t] == model.inf_b4treat[i, j, t] * (1 - 0.5*model.level1[i, j, t] )
model.treated_population = Constraint(model.Iset, model.Jset, model.Tset, rule=treatment_rule)
INFO: CONSTR treated_population[1,1,1]: -21.71460181746589 != 0.0
INFO: CONSTR treated_population[1,2,2]: -6.629094367820306 != 0.0
INFO: CONSTR treated_population[2,1,2]: -7.140628683877534 != 0.0
INFO: CONSTR treated_population[2,2,1]: -15.008019076531863 != 0.0

Output From Solver With NonNegativeReals boundaries applied + Code 413 Description from Knitro

Problem: 
- Lower bound: -inf
  Upper bound: inf
  Number of objectives: 1
  Number of constraints: 180
  Number of variables: 164
  Sense: unknown
Solver: 
- Status: warning
  Message: Knitro 11.1.0\x3a MIP\x3a All nodes have been explored. No integer feasible point found.; objective 1398.5939509578325; integrality gap 1.8e+308; 3 nodes; 104 subproblem solves
  Termination condition: maxIterations
  Id: 413
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

Code 413 Description: All nodes have been explored. No integer feasible point was found. The MIP optimality gap has not been reduced below the specified threshold, but there are no more nodes to explore in the branch and bound tree. If the problem is convex, this could occur if the gap tolerance is difficult to meet because of bad scaling or roundoff errors, or there was a failure at one or more of the subproblem nodes. This might also occur if the problem is nonconvex.

In theory, none of my variables should be negative when the model gets solutions. However, the question becomes the following: Why are the constraints becoming infeasible when I specify that all of my variables exist in the NonNegativeReals but are feasible when all domains are Real? Perhaps it has to do with the initial input values? Alternatively, perhaps there is a way that I can tell the model that if the model.susceptible variable goes negative, that instead it will just floor it to zero instead.

$\endgroup$
8
  • $\begingroup$ Would you mind adding an algebraic formulation of your model, or at least of the offending constraints? That might make it easier for folks here to follow your code. $\endgroup$ Oct 4, 2019 at 0:19
  • $\begingroup$ Will do. I will have to add it tomorrow because it will take some time to get the mathjax equations figured out! $\endgroup$ Oct 4, 2019 at 2:02
  • $\begingroup$ D.Gray Thanks for offering to add more equations written in MathJax, we generally expect people to do that (use MathJax) once they are around 500 reputation. IF they are particularly complicated, and it comes down to it, you could add the equations in image format and one of us will MathJaxify them for you. $\endgroup$
    – Rob
    Oct 4, 2019 at 11:33
  • 1
    $\begingroup$ Try following the advice in yalmip.github.io/debugginginfeasible Haven't really looked at your model, but why soukd it come as a shock that a model might be feasible, and then when additional constraints, such as variable nonnegativity are added, the model might become infeasible? Also, you state the specific parameter values used here are place holders - well, the feasibility or not of the model can be affected by the parameter values. On CVX forum, I've seen many times where unrealistic (random) made up values are used and the model is infeasible, but the original paper "worked" $\endgroup$ Oct 4, 2019 at 14:36
  • 1
    $\begingroup$ Not sure if this might be going on under the hood for KNITRO solving MINLP, but for continuous NLP applied to non-convex constraints, KNITRO can report local infeasiibility, which means it diidn't find a feasible solution from the starting point it used (was supplied); however, it might be able to find feasible solution from different starting point. Also, ability to find feasible solution depends not only on starting point, but on solver, algorithm choice in solver (KNTRO has several),& algorithm parameter values. In such case, reported infeasibility may not mean problem actually infeasible. $\endgroup$ Oct 4, 2019 at 15:38

1 Answer 1

6
$\begingroup$

There is a systematic way of finding the infeasibility of your problem. You would like to find the Irreducibly Inconsistent System (IIS) of your model. This technology is available in CPLEX and Gurobi for MIP and in BARON for MINLP.

Since you have implemented your model in Pyomo (in case you do not have a BARON license), you can submit the problem to the NEOS server for them to run the problem for you. You can either use the Pyomo-GAMS interface to generate a .gms file (make sure to use the symbolic_solver_labels=True, keepfiles=True option to keep the original names in the flattened GAMS model) as described here, or generate the .bar file to submit it to NEOS. If you have a BARON license run it directly.

Once you do that, you need to use the option CompIIS 1 for either GAMS or BARON directly to obtain the IIS for your problem. The formulation looks correct in my opinion, and from my own experience, the infeasibility usually comes from typos/data errors. I highly recommend the IIS computation to debug your models.

This procedure tends to be more expensive than the solution of the original MINLP, so consider using a small instance.

$\endgroup$
4
  • $\begingroup$ I'll have to check out what ISS is, but thank you. I can tell you where the infeasibility came from. The issue is with S_t; in that last term the negative variable is growing faster and if you observe the case where there is no non-negativity constraint, the 'S' variable becomes negative. When I implemented the non-negativity the model became infeasible because I was violating the non-negativity constraint. SO! I instead formulated the 'S_t' constraint as a Max{} and that solved the problem. $\endgroup$ Oct 10, 2019 at 18:16
  • 1
    $\begingroup$ @D.Gray It would be great if you could post your solution as an answer, to make it easier to find for future visitors. $\endgroup$ Oct 14, 2019 at 0:46
  • 1
    $\begingroup$ Ah yes, sorry, it slipped away from me. I'll get it posted for the a.m $\endgroup$ Oct 14, 2019 at 2:57
  • $\begingroup$ @D.Gray, did you have a chance to check IIS? I have a similar problem with my MINLP... $\endgroup$ Dec 3, 2019 at 17:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.