Simple problem with non-optimal solution

[Plug11]

I have a problem which I have simplified to a simple case with two continuous and one binary variable. There are no constraints (other than variable bounds) and while the objective is nonconvex, it is monotonic. I get the correct optimal solution with bonmin and MindtPy, but not consistently with SHOT. When I say not consistently, I mean that the results are sensitive to the exact algebraic representation I use in Pyomo.

I am very interested in using SHOT both for convex and nonconvex MINLP problems. Any thoughts on this would be greatly appreciated :-)

I am using a build of the following commit:

commit 033622c6f9989d48caf7eefe430d292286544711
Merge: 77da53f edd04d7
Author: Stefan Vigerske <[email protected]>
Date:   Mon Jun 27 14:50:10 2022 +0200

and I tried using both CBC and CPLEX MIP solvers (with slightly different results).

The following is the output from the different solvers (objective and variable values), for various different representations of the same problem. Note: when y==1.0 [optimum value], the value of the objective does not depend on x1. However this does not seem be the cause of this issue. I tried adding a small offset, so that the objective is always influenced by x1 and this caused bonmin and MindtPy to always give the same value for x1, as expected, but the behaviour of SHOT was still the same):

{'min': False, 'divide': False}
    SHOT_CBC:  OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: 6.5, x1: 0.6, x2: 0.2, y: 1.0
    bonmin:    OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': False, 'divide': True}
    SHOT_CBC:  OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: 6.5, x1: 0.6, x2: 0.2, y: 1.0
    bonmin:    OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
{'min': True, 'divide': False}
    SHOT_CBC:  OBJ: -2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: -2.0, x1: 0.5, x2: 0.1, y: 0.0
    bonmin:    OBJ: -13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': True, 'divide': True}
    SHOT_CBC:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    bonmin:    OBJ: -13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: -13.0, x1: 0.6, x2: 0.1, y: 1.0

and this is the self-contained code to replicate:

import pyomo.environ as pyo
import pyomo.contrib.mindtpy.MindtPy as mindtpy


def create_model(min=False, divide=False):

    model = pyo.ConcreteModel()

    model.x1 = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 0.6))
    model.x2 = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.1, 0.2))
    model.y = pyo.Var(domain=pyo.Boolean, initialize=1)

    y_mult = 1.3
    if divide:
        obj = (1 - model.y) / model.x1 + model.y * y_mult / model.x2
    else:
        obj = (1 - model.y) * (model.x1 ** -1) + model.y * y_mult * (model.x2 ** -1)

    if min:
        model.OBJ = pyo.Objective(expr = -1 * obj, sense=pyo.minimize)
    else:
        model.OBJ = pyo.Objective(expr = obj, sense=pyo.maximize)

    return model


solvers = {}

shotCbc = pyo.SolverFactory("SHOT")
shotCbc.options["Dual.MIP.Solver"] = 2
solvers["SHOT_CBC"] = shotCbc

shotCplex = pyo.SolverFactory("SHOT")
shotCplex.options["Dual.MIP.Solver"] = 0
solvers["SHOT_CPLEX"] = shotCplex

solvers["bonmin"] = pyo.SolverFactory("bonmin")

mindtpy_solver = mindtpy.MindtPySolver()
mindtpy_solver.CONFIG["mip_solver"] = "cbc"
solvers["MindtPy"] = mindtpy_solver

params = [
    {"min": False, "divide": False}, 
    {"min": False, "divide": True}, 
    {"min": True, "divide": False}, 
    {"min": True, "divide": True}, 
]
results = {}
for p in params:
    cur_res = {}
    for name, sol in solvers.items():
        model = create_model(**p)
        sol.solve(model, tee=True)
        cur_res[name] = model
    results[str(p)] = cur_res

for p, res in results.items():
    print(p)
    for name, r in res.items():
        print(f"    {name}:{' '*(10-len(name))}OBJ: {pyo.value(r.OBJ):.1f}, x1: {pyo.value(r.x1):.1f}, x2: {pyo.value(r.x2):.1f}, y: {pyo.value(r.y):.1f}")

[andreaslundell]

Hi,

Unfortunately I do not currently have SHOT setup with Pyomo, so I cannot test your problem easily. But do you think you could rerun SHOT with Output.Debug.Enable=true? This will create a folder with debug information (e.g. the subproblems). The path to this folder can be set with Output.Debug.Path=/path/to/folder. If you can attach these files (as a zip-file) to the issue I could take a look. You can also send them by e-mail to me at [email protected].

Regards, Andreas

[grahamsparrow-8451]

Hi Andreas,

Thanks for replying and offering to have a look at this. I am attaching the log files. There are 4 separate runs corresponding to the different representations in the pyomo code above and in subdirectories with matching names. Please let me know if there is anything else that would be useful. I am happy to run further tests and also to look at the C++ code, if this could help; I have experience of writing numerical C++ code, but am not familiar with the SHOT codebase.

SHOT_issue166_logs.zip

Best regards,
Graham

[imliuyifan]

I tried to replicate with Graham's code with :

python: 3.11.4
Pyomo version: 6.6.1

Shot is using nlp solver Ipopt: 3.14.10 lp solver glpk: 5.0
Bonmin 1.8.9 is using Cbc 2.10.8 and Ipopt 3.12.13
and got these:

{'min': False, 'divide': False}
    SHOT_CBC:  OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: 6.5, x1: 0.6, x2: 0.2, y: 1.0
    bonmin:    OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': False, 'divide': True}
    SHOT_CBC:  OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: 6.5, x1: 0.6, x2: 0.2, y: 1.0
    bonmin:    OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': True, 'divide': False}
    SHOT_CBC:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    bonmin:    OBJ: -13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': True, 'divide': True}
    SHOT_CBC:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    bonmin:    OBJ: -13.0, x1: 0.6, x2: 0.1, y: 1.0
    MindtPy:   OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0

[andreaslundell]

This is indeed a bug in SHOT when maximizing. I have not located it yet, but will when I return to work in a couple of weeks. In the meantime, I suggest minimizing :-)

[Plug11]

Thanks for the update, @andreaslundell. It is interesting that you have found that the issue relates to maximizing; @imliuyifan also found that he only got the incorrect results when maximizing. However, in my tests, one of the minimizing cases also gave the wrong answer. I will look into my code/results again.

[andreaslundell]

Sorry for the delay, but I was finally able to sit down and locate the bug.

Since this is a nonconvex problem, SHOT does not guarantee to find the optimal solution (13 or -13). It will return the best solution found, and in the case of SHOT's AMPL interface, which Pyomo uses, it will return a message Solved to local optimality if the solution can not be proven to be global (dual bound = primal bound). However, there was a bug in SHOT when maximizing a nonlinear objective function, that prevented it from working as expected.

This should now be fixed in the branch https://github.com/coin-or/SHOT/tree/issue166. I modified your example somewhat to also test with Gurobi. I got the following output:

{'min': False, 'divide': False}
    SHOT_GUR:  OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
    SHOT_GUR_2:OBJ: 13.0, x1: 0.6, x2: 0.1, y: 1.0
    SHOT_CBC:  OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CBC_2:OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: 6.5, x1: 0.6, x2: 0.2, y: 1.0
    SHOT_CPLEX_2:OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': False, 'divide': True}
    SHOT_GUR:  OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_GUR_2:OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CBC:  OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CBC_2:OBJ: 2.0, x1: 0.5, x2: 0.1, y: 0.0
    SHOT_CPLEX:OBJ: 6.5, x1: 0.6, x2: 0.2, y: 1.0
    SHOT_CPLEX_2:OBJ: 13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': True, 'divide': False}
    SHOT_GUR:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_GUR_2:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CBC:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CBC_2:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX_2:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
{'min': True, 'divide': True}
    SHOT_GUR:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_GUR_2:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CBC:  OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CBC_2:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0
    SHOT_CPLEX_2:OBJ: -13.0, x1: 0.5, x2: 0.1, y: 1.0

So, as we can see, SHOT with Gurobi manages to find the optimal solution also when maximizing. Here the SHOT_GUR_2, SHOT_CBC_2 and SHOT_CPLEX_2 runs are using SHOT's own functionality for solving NLP calls (Primal.FixedInteger.Solver=2) instead of Ipopt. This finds the solution (y=1) also for Cbc. Normally using the internal NLP call functionality in SHOT is not recommended for nonconvex problems, but in this case, it seems to work well.

Could you verify @Plug11 and @grahamsparrow-8451 that this works for you.