# Understanding why this MINLP formulation is infeasible

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
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]


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.

• 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. Oct 4 '19 at 0:19
• Will do. I will have to add it tomorrow because it will take some time to get the mathjax equations figured out! Oct 4 '19 at 2:02
• 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.
– Rob
Oct 4 '19 at 11:33
• 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" Oct 4 '19 at 14:36
• 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. Oct 4 '19 at 15:38

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.