Integer Decision Variables Always Forced to Zero in Minimization Problem (MINLP)

Question

I am trying to write out a MINLP problem of optimal control for an invasive species and the code that I have for my PYOMO model is below. Some of the initialization values take from an Excel spreadsheet, so they are not made explicit in the formulation below. Additionally, these values aren't real data, but simulated data I have made just to try and get my model running.

The problem that I am encountering is that my binary variable model.level1, which emphasizes whether or not a cell receives treatment, is always slammed to zero in my objective function which indicates that no control is applied to any cells and this effect cascades throughout the rest of the equations and results in unrealistic results even for fake data.

Perhaps I am overlooking something in my constraints that are causing the binary variable to be slammed to 0 rather than having a variety of different 0's and 1's throughout the solution?

model = ConcreteModel()
Imax = 2
Jmax = 2
Tmax = 3
Kmax = 2
model.Iset = RangeSet(1, Imax)  # e.g. i = {1, 2, 3}
model.Jset = RangeSet(1, Jmax)
model.Tset = RangeSet(1, Tmax)
model.Kset = RangeSet(1, Kmax)

model.juvsurv = Param(initialize=0.60)
model.juvdeath = Param(initialize=0.40)
model.lambd = Param(initialize=0.20)
model.alpha = Param(initialize=2)

model.juvenille = Var(model.Iset, model.Jset, model.Tset, model.Kset, within=NonNegativeReals, initialize=initial_values_ext(3, Tmax, 2, juvenille_init))
model.susceptible = Var(model.Iset, model.Jset, model.Tset,within=NonNegativeReals, initialize=initial_values(3, Tmax, susceptible_init))
model.juvenilleTotal = Var(model.Iset, model.Jset, model.Tset, within=NonNegativeReals)
model.inf_b4treat = Var(model.Iset, model.Jset, model.Tset, within=NonNegativeReals, initialize=initial_values(3, Tmax, inf_b4treat_init))
model.inf_treated = Var(model.Iset, model.Jset, model.Tset, within=NonNegativeReals)
model.obj = Var(model.Iset, model.Jset, model.Tset, within=NonNegativeReals)

def objective_rule(model):
    return sum (sum (sum (model.obj[i,j,t] for i in model.Iset ) for j in model.Jset ) for t in model.Tset )
model.damages = Objective(rule=objective_rule, sense=minimize)

def obj_rule(model, i,j,t):
    return model.obj[i,j,t] == 10*model.inf_b4treat[i,j,t] + 5*model.level1[i,j,t]*model.inf_b4treat[i,j,t]
model.object = Constraint(model.Iset, model.Jset, model.Tset, rule=obj_rule)

# Constraint 1: juvenilles that advance to the next age class (eq. 1)
def juv_advance_rule(model, i, j, t, k):
    if t != Tmax and k != Kmax:
        return model.juvenille[i, j, t + 1, k + 1] == model.juvenille[i, j, t, k] * model.juvsurv
    return Constraint.Skip
model.juv_advance = Constraint(model.Iset, model.Jset, model.Tset, model.Kset, rule=juv_advance_rule)

# Constraint 2: total number of juvenilles in all age classes on cell (i,j)
def juv_total_rule(model, i, j, t):
    return model.juvenilleTotal[i, j, t] == sum(model.juvenille[i, j, t, k] for k in model.Kset)
model.juv_total = Constraint(model.Iset, model.Jset, model.Tset, rule=juv_total_rule)

# Constraint 3: recruitment of seedlings to the first juvenille age class.
def juv_recruit_rule(model, i, j, t, k):
    if k == 1 and t != Tmax:
        return model.juvenille[i, j, t + 1, k] == model.juvdeath * model.juvenilleTotal[i, j, t]
    else:
        return Constraint.Skip
model.juv_recruit = Constraint(model.Iset, model.Jset, model.Tset, model.Kset, rule=juv_recruit_rule)

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

# Constraint 5(10): Population Growth:
def infested_growth_rule(model, i, j, t):
    if t == Tmax:
        return Constraint.Skip
    else:
        return model.inf_b4treat[i, j, t + 1] == ( model.inf_treated[i, j, t]**2 / (model.inf_treated[i,j,t]**2 + model.alpha) ) * model.susceptible[i,j,t] + 10
model.inf_growth = Constraint(model.Iset, model.Jset, model.Tset, rule=infested_growth_rule)

# Constraint 9: Treated Infestation
def treatment_rule(model, i, j, t):
    return model.inf_treated[i, j, t] == model.inf_b4treat[i, j, t] * (1 - model.level1[i, j, t] )
model.treated_pop = Constraint(model.Iset, model.Jset, model.Tset, rule=treatment_rule)

def budget_rule(model):
    return sum(sum(sum(2*model.level1[i,j,t] for i in model.Iset) for j in model.Jset) for t in model.Tset) <= 3
model.budget = Constraint(model.Iset, model.Jset, model.Tset, rule=budget_rule)

```

Answer 1

You can try adding a constraint forcing one of the affected variable to be nonzero. If the model becomes infeasible, you can try to find the conflicting constraints. If the model stays feasible, this means that your objective function represents other priorities than you expected.

Answer 2

I only skimmed your model so others may be better able to point to the error directly, but here are some reasons this may occur:

• Constraints: as you mention, perhaps they're set so that it's not possible to apply the treatment. E.g, is the budget accidentally set too tightly so it can't afford it? (There are several sums in your budget_rule - I'd double check that constraint)
• Objective function: perhaps it hurts the OF value when treatment is applied. E.g., does the cost of treatment outweigh the benefit? Also I'm not sure why model.inf_b4treat is included in both the not treated and treated terms (could be fine). Also, always good to double-check whether you mean to minimize vs. maximize.
• Inputs: is the problem initialized so that no one needs treatment?

Other ideas to debug:

• Test the model with the an extreme inputs and/or adjusted/removed constraints. E.g., no budget limit, no treatment cost
• If still not clear, iteratively remove constraints and resolve. See what the limiting one may be.

Best wishes!