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?

Use a global MINLP solver instead. – Erwin Kalvelagen Apr 14 '23 at 05:23 — Erwin Kalvelagen, Apr 14 '23 at 05:23

AirSquid · Answer 1 · 2023-04-09T16:30:50.990

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

You still had multiplied variables in constraint 1 for no reason.
No need to initialize variables on this small of a model, esp. when the initialization is invalid
print the model and review it, same w/ result as shown
If this gets bigger, I'd index P and y instead of hand-jamming them all individually.
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))

I've made the changes - removing con2 - con6 constraints and removing y from the objective function. However, it's now converging on a different solution, choosing P2 and P4, strangely. Setting tee=True, it's saying that this is Optimal. Here is the colab: https://colab.research.google.com/drive/1UuIXdnVjQTEMyWKzS5P4bu10K1uuVPFQ?usp=sharing — wizardjoe, Apr 09 '23 at 15:50
Also, if we're making it linear, then it will be similar to the other LP problem you helped me solve (in PuLP). However I don't know how to change that other problem to make it work with P in the objective function instead of y, where P is multiplied by the weights: https://stackoverflow.com/questions/75955386/pulp-error-when-constraining-variables-that-are-mutually-exclusive/75959614?noredirect=1 — wizardjoe, Apr 09 '23 at 15:52

Marcin Orlowski · Answer 2 · 2023-04-09T11:11:03.977

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.

I made the update and ran it, but now the solver has set all ```yi``` variables to 0 — wizardjoe, Apr 09 '23 at 07:43
It seems that the solver is getting stuck in a local optimum or failing to find a feasible solution. See answer edit. — Marcin Orlowski, Apr 09 '23 at 08:08
The new constraint didn't work. But I modified the problem and got a little further. I replaced the bounds with constraints, and also changed the bounds to make the problem more apparent. With the new arrangement, I don't need to force the sum of the ```yi``` variables to be >= 1. However, it's still arriving at a suboptimal solution. This time the correct answer should be y5 = 1 and p5 = 50, but it's doing y3 and y4, with P3 = 29.999 and P4 = 20. — wizardjoe, Apr 09 '23 at 08:18

Pyomo MINLP solver not picking the optimal outcome

2 Answers2

EDIT: code added.