cvxpy: Linear mixed-integer solver problem is failing with GLPK-MI

Update: see also comment below.

I’m having similar issues with a similar problem, I can also reproduce the problem in the bug description. After a bit of trial and error, I suspect that when using GLPK_MI, equality constraints are not supported or not working correctly. There might be an easy workaround by adding the constraint with inequalities in both directions. At least for your example it seems to spit out the same result as it does with Gurobi. Disclaimer: lots of guess work involved here 😃

Environment

>>> cvxpy.__version__
'1.1.3'
>>> cvxopt.__version__

glpk=4.65=he80fd80_1002

$ /opt/gurobi810/linux64/bin/gurobi_cl --version
Gurobi Optimizer version 8.1.0 build v8.1.0rc1 (linux64)

Complete cona environment

Minimal example 1

the original example with making the solver choice explicit via solver=cp.GLPK_MI.

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.],
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) == 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve(solver=cp.GLPK_MI, verbose=True)
print(f"Optimal disorder is {optimal_disorder}")

Output:

Traceback (most recent call last):
  File "/home/p/abcvoting/survey/example05-pav.py", line 16, in <module>
    optimal_disorder = prob.solve(solver=cp.GLPK_MI, verbose=True)
  File "/opt/miniconda3/envs/abcvoting/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 395, in solve
    return solve_func(self, *args, **kwargs)
  File "/opt/miniconda3/envs/abcvoting/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 747, in _solve
    self.unpack_results(solution, solving_chain, inverse_data)
  File "/opt/miniconda3/envs/abcvoting/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 1055, in unpack_results
    "Try another solver, or solve with verbose=True for more "
cvxpy.error.SolverError: Solver 'GLPK_MI' failed. Try another solver, or solve with verbose=True for more information.

Process finished with exit code 1

Minimal example 2

same as example 1, but replacing the constraint with two inequalities (and adding more output).

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.],
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) <= 1, cp.matmul(A, x) >= 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve(solver=cp.GLPK_MI, verbose=True)

print(f"Optimal disorder is {optimal_disorder}")
print(f"x: {x.value}")

Output:

/opt/miniconda3/envs/abcvoting/bin/python /home/p/abcvoting/survey/example05-pav.py
      0: obj =   0.000000000e+00 inf =   5.000e+00 (5)
      5: obj =   3.360000000e+00 inf =   0.000e+00 (0)
*     9: obj =   3.360000000e+00 inf =   0.000e+00 (0)
+     9: mip =     not found yet >=              -inf        (1; 0)
+     9: >>>>>   3.360000000e+00 >=   3.360000000e+00   0.0% (1; 0)
+     9: mip =   3.360000000e+00 >=     tree is empty   0.0% (0; 1)
Optimal disorder is 3.36
x: [1. 0. 0. 0. 0. 0. 1. 1. 1. 0.]

Process finished with exit code 0

Minimal example 3

Same thing with Gurobi and equality constraint for cross-checking the solution:

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.],
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) == 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve(solver=cp.GUROBI, verbose=True)

print(f"Optimal disorder is {optimal_disorder}")
print(f"x: {x.value}")

Output:

Parameter OutputFlag unchanged
   Value: 1  Min: 0  Max: 1  Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0  Min: 0  Max: 1  Default: 0
Optimize a model with 5 rows, 10 columns and 13 nonzeros
Variable types: 0 continuous, 10 integer (10 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [4e-01, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 5.0000000
Presolve removed 5 rows and 10 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 8 available processors)

Solution count 2: 3.36 5 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.360000000000e+00, best bound 3.360000000000e+00, gap 0.0000%
Optimal disorder is 3.36
x: [1. 0. 0. 0. 0. 0. 1. 1. 1. 0.]

Process finished with exit code 0

We’re pushing SCIP now: https://www.cvxpy.org/install/#install-with-scip-support It’s unfortunate that people now need to install a solver to solve mixed integer problems. The issue is that the solutions from ECOS_BB were often totally wrong.

It’s actually pretty simple: cvxopt.glpk.ilp() doesn’t support linear programs without inequality constraint, so it’s not really a bug. CVXPY could probably pass on the error message to the user in a nicer format, but I wouldn’t say this is a bug. But it should be possible to add support for it easily, I’m just not sure what’s the right way to go here. Is there any reason, why one would require at least one in-equality constraint for the solver?

I guess we could simply use the equality constraints also for the in-equality constraints if there is no in-equality constraint, i.e. add the workaround from my previous post to CVXPY. Or we could add just one (arbitrary?) in-equality constraint as equality constraint. I guess solving a linear program without any constraints, doesn’t make sense, right? So there must be some constraint. Is there a performance or numerical difference which constraint is chosen? Is there any better solution than this? Maybe it needs to be fixed upstream in CVXOPT?

It might be helpful to improve error handling es well. Debugging might be easier, if a an error such as “m must be a positive integer” is passed to the user.

The docstring of cvxopt.glpk.ilp():

Solves a mixed integer linear program using GLPK.

(status, x) = ilp(c, G, h, A, b, I, B)

PURPOSE
Solves the mixed integer linear programming problem

    minimize    c'*x
    subject to  G*x <= h
                A*x = b
                x[k] is integer for k in I
                x[k] is binary for k in B

ARGUMENTS
c            nx1 dense 'd' matrix with n>=1

G            mxn dense or sparse 'd' matrix with m>=1

h            mx1 dense 'd' matrix

A            pxn dense or sparse 'd' matrix with p>=0

b            px1 dense 'd' matrix

I            set of indices of integer variables

B            set of indices of binary variables

I’m new to the internals of CVXPY, but would be happy to provide a patch if somebody could help me a bit along the way.

Woah great work @lumbric ! That’s actually quite an hacky yet elegant way to make this work. I’ll try to see if my unit tests fail or not when I’m using this workaround.