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)
```