Pyomo MINLP solver not picking the optimal outcome

Question

I'm using Pyomo to solve an MINLP problem. It seems to run without issue, but the solution it's picking does not seem to be optimal.

The gist of the problem is that I'm trying to optimize a function of the form (C1 * P1 * y1) + (C2 * P2 * y2)...(C5 * P5 * y5)= Z, where C1...C5 are constants, and P1...P5 and y1...y5 are decision variables. The P variables are continuous and the y variables are binary, with the y variables functioning as selectors for the P variables (ie, if y5 is 1, then P5 is factored into the objective, etc). We are also given a constraint that the sum of P1y1 + ... P5y5 <= 50.

Here's the code:

from pyomo.environ import *

# create a ConcreteModel object
model = ConcreteModel()

# define the decision variables
model.P1 = Var(within=NonNegativeReals, initialize=5.0)
model.P2 = Var(within=NonNegativeReals, initialize=15.0)
model.P3 = Var(within=NonNegativeReals, initialize=25.0)
model.P4 = Var(within=NonNegativeReals, initialize=35.0)
model.P5 = Var(within=NonNegativeReals, initialize=7.5)
model.y1 = Var(within=Binary, initialize=0)
model.y2 = Var(within=Binary, initialize=0)
model.y3 = Var(within=Binary, initialize=0)
model.y4 = Var(within=Binary, initialize=0)
model.y5 = Var(within=Binary, initialize=0)

# define the objective function
model.obj = Objective(expr=1*model.P1*model.y1 + 1.5*model.P2*model.y2 + 1.7*model.P3*model.y3 +
                      2*model.P4*model.y4 + 2.5*model.P5*model.y5, sense=maximize)

# define the constraints
model.con1 = Constraint(expr=model.P1*model.y1 + model.P2*model.y2 + model.P3*model.y3 +
                        model.P4*model.y4 + model.P5*model.y5 <= 50)
model.con2 = Constraint(expr=model.P1 <= 1000000 * model.y1)
model.con3 = Constraint(expr=model.P2 <= 1000000 * model.y2)
model.con4 = Constraint(expr=model.P3 <= 1000000 * model.y3)
model.con5 = Constraint(expr=model.P4 <= 1000000 * model.y4)
model.con6 = Constraint(expr=model.P5 <= 1000000 * model.y5)
model.con7 = Constraint(expr=model.y1 + model.y2 + model.y3 <= 1)
model.con8 = Constraint(expr=model.y4 + model.y5 <= 1)

# add upper and lower bounds for P variables
model.con9 = Constraint(expr=model.P1 >= model.y1 * 10)
model.con10 = Constraint(expr=model.P1 <= model.y1 * 20)

model.con11 = Constraint(expr=model.P2 >= model.y2 * 20)
model.con12 = Constraint(expr=model.P2 <= model.y2 * 30)

model.con13 = Constraint(expr=model.P3 >= model.y3 * 30)
model.con14 = Constraint(expr=model.P3 <= model.y3 * 500)

model.con15 = Constraint(expr=model.P4 >= model.y4 * 15)
model.con16 = Constraint(expr=model.P4 <= model.y4 * 30)

model.con17 = Constraint(expr=model.P5 >= model.y5 * 30)
model.con18 = Constraint(expr=model.P5 <= model.y5 * 500)

# specify the solver to use
solver = SolverFactory('mindtpy') 

# solve the problem
solver.solve(model, strategy="OA", mip_solver='glpk', nlp_solver='ipopt')

# print the results
print('P1 =', value(model.P1))
print('P2 =', value(model.P2))
print('P3 =', value(model.P3))
print('P4 =', value(model.P4))
print('P5 =', value(model.P5))
print('y1 =', value(model.y1))
print('y2 =', value(model.y2))
print('y3 =', value(model.y3))
print('y4 =', value(model.y4))
print('y5 =', value(model.y5))

And here's the result:

P1 = 9.727046454907866e-13
P2 = 9.727046454907866e-13
P3 = 29.999999708353013
P4 = 20.000000790394036
P5 = 9.727046454907866e-13
y1 = 0.0
y2 = 0.0
y3 = 1.0
y4 = 1.0
y5 = 0.0

Note that the solver has switched on y3 and y4, and given P3 and P4 values. If we do the math, the objective function will equal around 91. However, the optimal solution should be y5 = 1 and P5 = 50, to yield an objective function value of 125.

What's wrong here?

Answer 1

The constraints con2, con3, con4, con5 and con6 don't work as intended. Try replacing these constraints with:

model.con2 = Constraint(expr=model.P1 <= 1000000 * model.y1)
model.con3 = Constraint(expr=model.P2 <= 1000000 * model.y2)
model.con4 = Constraint(expr=model.P3 <= 1000000 * model.y3)
model.con5 = Constraint(expr=model.P4 <= 1000000 * model.y4)
model.con6 = Constraint(expr=model.P5 <= 1000000 * model.y5)

This ensures that the continuous variables (P1, P2, P3, P4, P5) are effectively "switched off" when the corresponding binary variables (y1, y2, y3, y4, y5) are equal to zero. In your original constraints (con2 to con6), you multiplied the continuous variables (Pi) by the binary variables (yi), which doesn't effectively enforce the desired behavior. The original constraints only ensure that if yi == 1, then Pi can take any value within its bounds, but they don't enforce Pi to be zero when yi == 0.

---

EDIT #1: You can add an additional constraint to your model to help guide the solver towards a better solution. You mentioned that the optimal solution should have y4 == 1 and P4 == 50, so we can set an additional constraint that forces at least one of the y variables to be equal to 1:

model.con11 = Constraint(expr=model.y1 + model.y2 + model.y3 + model.y4 + model.y5 >= 1)

---

EDIT #2: You can increase the maximum number of iterations for the solver, which would give it more time to find a better solution:

solver.solve(model, strategy='OA', mip_solver='glpk', nlp_solver='ipopt', max_iter=5000)

Another option worth considering would be creating an additional constraint to prevent the current sub-optimal solution. For example, if the current solution is y3 == 1 and y4 == 1 with P3 == 29.999 and P4 == 20, you can add a constraint that prohibits this combination and force solver to search for other solutions:

model.con12 = Constraint(expr=model.y3 + model.y4 <= 1)

Finally, if still unsuccessful, maybe see if other solvers like bonmin, couenne or baron would be more handy.

Answer 2

As your model is written, after the edits, constraints 2-6 are superfluous. The other later constraints clamp the P variables to zero when y is zero.

So the P variables are semi-continuous in this model, meaning they are either zero or within a range. If this is your desire, the model can simply be made linear by removing y variables from your objective function as they are also superfluous. If y is zero, so is the associated P leaving you with an objective function that is the sum(P * weight).

Lastly, I don’t see in your code where you are inspecting the solution object to ensure that the status is “OPTIMAL”. That is a huge mistake. The results are gibberish unless the solver claims optimal.

EDIT: code added.

from the last iteration...

1. You still had multiplied variables in constraint 1 for no reason.

2. No need to initialize variables on this small of a model, esp. when the initialization is invalid

3. print the model and review it, same w/ result as shown

4. If this gets bigger, I'd index P and y instead of hand-jamming them all individually.

5. If you are making a linear model, use a linear solver...if it barfs, something is wrong.

---

from pyomo.environ import *

# create a ConcreteModel object
model = ConcreteModel()

# define the decision variables
model.P1 = Var(within=NonNegativeReals)
model.P2 = Var(within=NonNegativeReals)
model.P3 = Var(within=NonNegativeReals)
model.P4 = Var(within=NonNegativeReals)
model.P5 = Var(within=NonNegativeReals)
model.y1 = Var(within=Binary)
model.y2 = Var(within=Binary)
model.y3 = Var(within=Binary)
model.y4 = Var(within=Binary)
model.y5 = Var(within=Binary)

# define the objective function
model.obj = Objective(expr=1*model.P1 + 1.5*model.P2 + 1.7*model.P3 +
                      2*model.P4 + 2.5*model.P5, sense=maximize)

# define the constraints
model.con1 = Constraint(expr=model.P1 + model.P2 + model.P3 +
                        model.P4 + model.P5 <= 50)

model.con7 = Constraint(expr=model.y1 + model.y2 + model.y3 <= 1)
model.con8 = Constraint(expr=model.y4 + model.y5 <= 1)

# add upper and lower bounds for P variables
model.con9 = Constraint(expr=model.P1 >= model.y1 * 10)
model.con10 = Constraint(expr=model.P1 <= model.y1 * 20)

model.con11 = Constraint(expr=model.P2 >= model.y2 * 20)
model.con12 = Constraint(expr=model.P2 <= model.y2 * 30)

model.con13 = Constraint(expr=model.P3 >= model.y3 * 30)
model.con14 = Constraint(expr=model.P3 <= model.y3 * 500)

model.con15 = Constraint(expr=model.P4 >= model.y4 * 15)
model.con16 = Constraint(expr=model.P4 <= model.y4 * 30)

model.con17 = Constraint(expr=model.P5 >= model.y5 * 30)
model.con18 = Constraint(expr=model.P5 <= model.y5 * 500)


# QA the model:
model.pprint()

# specify the solver to use
solver = SolverFactory('cbc') 

# solve the problem
result = solver.solve(model) # strategy="OA", mip_solver='glpk', nlp_solver='ipopt', tee=True, nlp_solver_tee=True, mip_solver_tee=True)

# review the solution and solve status:
print(result)


# print the results
print('P1 =', value(model.P1))
print('P2 =', value(model.P2))
print('P3 =', value(model.P3))
print('P4 =', value(model.P4))
print('P5 =', value(model.P5))
print('y1 =', value(model.y1))
print('y2 =', value(model.y2))
print('y3 =', value(model.y3))
print('y4 =', value(model.y4))
print('y5 =', value(model.y5))