9
$\begingroup$

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)

```
$\endgroup$
1
  • 1
    $\begingroup$ It might help to write a minimal, reproducible example (reprex). This will make it easier for folks here to diagnose the problem, and a bonus is that the process of building the reprex is in some ways similar to the process of debugging, so you might actually find the bug on your own. $\endgroup$ Sep 12, 2019 at 21:31

2 Answers 2

5
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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!

$\endgroup$
1
  • $\begingroup$ Initially yes I set the problem up so that model.level1 is 0. My thought process here was that running the model would provide the optimal values for that variable. It seemed strange to me that some areas would have the model.level1 variable be 1 while others would be 0 for no particular reason. I will give these debugging methods a go and see what I can figure out. $\endgroup$ Sep 11, 2019 at 17:48

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.