Optimal solution found, but constraint is violated

[kuvychko]

Hello!

I found this issue while playing with a toy problem (finding greatest common divisor using MIP). This formulation works fine with GLPK (optimal solution is found and it is correct). When I use CBC the solver flies off and reports an incorrect solution that violates a constraint but termination message says "Model was solved to optimality (subject to tolerances), and an optimal solution is available." If I add additional constraints on decision variables (to prevent "fly-off"), CBC converges on a correct solution. The problem persists across the range of inputs for GCD (CBC fails for all inputs I tried, while GLPK finds correct optimal solutions). Thank you for your help with this!!

Here is a control file (generated using Pyomo, looking for GCD(9, 27). CBC gives a solution of -900:

'a': 9,
'b': 27,
'A': -12345679000,
'B': 4115226300,
'gcd': -900.0,

2 Var Declarations
    A : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :  None :  None :  None : False :  True : Integers
    B : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :  None :  None :  None : False :  True : Integers

1 Objective Declarations
    objective : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 9*A + 27*B

1 Constraint Declarations
    constraint : Size=1, Index=None, Active=True
        Key  : Lower : Body       : Upper : Active
        None :   1.0 : 9*A + 27*B :  9.0 :   True

4 Declarations: A B objective constraint

If I add extra constraints on decision variables to prevent "fly-off", CBC finds the right solution. Here is the corresponding control file (same problem GCD(9,27)):

'a': 9,
'b': 27,
'A': -26.0,
'B': 9.0,
'gcd': 9.0,

2 Var Declarations
    A : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :   -27 :  None :    27 : False :  True : Integers
    B : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :    -9 :  None :     9 : False :  True : Integers

1 Objective Declarations
    objective : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 9*A + 27*B

1 Constraint Declarations
    constraint : Size=1, Index=None, Active=True
        Key  : Lower : Body       : Upper : Active
        None :   1.0 : 9*A + 27*B :  9.0 :   True

4 Declarations: A B objective constraint

Pyomo model:

import time
import pyomo.environ as pyo

def gcd_ip1(a, b, solver='cbc', additional_constraints=False, filename=None):
    """Compute the greatest common divisor of a and b using MIP (alternative formulation)
    
    Uses the following model (alternative formulation using inequality constraint):
    min A*a + B*b
    s.t. A*a + B*b >= 1

    gcd(a, b) = A*a + B*b    
    """
    # Initialize the model
    model = pyo.ConcreteModel()

    # Define decision variables
    if additional_constraints:
        model.A = pyo.Var(domain=pyo.Integers, bounds=(-b, b))
        model.B = pyo.Var(domain=pyo.Integers, bounds=(-a, a))
    else:
        model.A = pyo.Var(domain=pyo.Integers)
        model.B = pyo.Var(domain=pyo.Integers)
        
    # Define the objective
    model.objective = pyo.Objective(
        expr=model.A*a + model.B*b, sense=pyo.minimize
    )

    # Define constraint
    model.constraint = pyo.Constraint(expr= (1, model.A*a + model.B*b, min(a, b)))

    # Write the model to a file
    if filename:
        with open(filename, 'w') as f:
            model.pprint(ostream=f)
    
    # Solve the model and time it
    solver = pyo.SolverFactory(solver)
    start_time = time.time()
    result = solver.solve(model)
    elapsed_time = time.time() - start_time

    # Compile the results into a dictionary and return
    results = {
        'a': a,
        'b': b,
        'A': pyo.value(model.A),
        'B': pyo.value(model.B),
        'gcd': pyo.value(model.A*a + model.B*b),
        'time': elapsed_time,
        'Status': result.solver.status,
        'Termination Condition': result.solver.termination_condition,
        'Termination Message': result.solver.termination_message,
        'Optimal Value (z)': model.objective()
    }

    return results

My version of CBC comes from Cbc-master-win64-msvc16-mt.zip which I got from https://www.coin-or.org/download/binary/Cbc/ It looks like a development version:

(base) C:\Users\igork>cbc --version
Welcome to the CBC MILP Solver
Version: devel
Build Date: Apr 27 2021

I'm using Pyomo 6.8.2.

[jjhforrest]

Cbc gets the correct solution - but does not print it correctly - in debug mode -
p primalColumnSolution[0]@2
$3 = {12345678898, -4115226299}
which is a valid answer. The trouble is that the format used is %15.8g.

Will change code to give better printing with very large integer variables

[kuvychko]

@jjhforrest - I see, thank you for looking into this!