5

I am using CVXOPT to solve a very simple problem:

min -7890424934354.171875*x1 -7890424934354.274414*x2 -7890424934354.246093*x3
s.t: 
  x1 + x2 + x3 = 1
  x1,x2,x3 are binary

We can see that the optimal solution should be obviously:

x1 =0; x2 = 1; x3 = 0

However I didn't get a correct answer using ILP from CVXOPT(I know the above problem is too simple to use ILP, but I am just curious). A detailed description about ILP of CVXOPT is here.

My program is like this:

from cvxopt.glpk import ilp
from cvxopt import matrix
c = matrix([-7890424934354.171875,-7890424934354.274414,-7890424934354.246093],tc='d')
G = matrix(0.0, (1,3)) #since I do not have a constraint like G*x <= h, I make them zeros here
h = matrix(0.0, (1,1))
A = matrix([1,1,1],tc='d')
b = matrix(1,tc='d')
(status, x) = ilp(c,G,h,A.T,b,B=set([0,1,2]))

The result is:

GLPK Integer Optimizer, v4.61
2 rows, 3 columns, 3 non-zeros
3 integer variables, all of which are binary
Preprocessing...
1 row, 3 columns, 3 non-zeros
3 integer variables, all of which are binary
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 1
Solving LP relaxation...
GLPK Simplex Optimizer, v4.61
1 row, 3 columns, 3 non-zeros
*     0: obj =  -7.890424934e+12 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Integer optimization begins...
+     0: mip =     not found yet >=              -inf        (1; 0)
+     0: >>>>>  -7.890424934e+12 >=  -7.890424934e+12   0.0% (1; 0)
+     0: mip =  -7.890424934e+12 >=     tree is empty   0.0% (0; 1)
INTEGER OPTIMAL SOLUTION FOUND
optimal
[ 0.00e+00]
[ 0.00e+00]
[ 1.00e+00]

which doesn't return the optimal solution.

But if I change my objective function to -171875*x1 - 274414*x2 - 246093 * x3, I can get a correct answer which is x1 = 0, x2 = 1, x3 = 0.

I am really confused why this happen:

I guessed firstly whether floating-point values like -7890424934354.171875 lose precision when passed to ILP, but it seems this is not the reason.

Another guess is that I shouldn't make G and h zeros if I don't have a constraint like G*x <= h. But what should I do if I want to use ILP when G and h are empty since it requires

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

I'd appreciate any advice. Many thanks in advance.

miguelmorin
  • 5,025
  • 4
  • 29
  • 64
Sue
  • 113
  • 1
  • 9
  • 2
    Coefficients should have reasonable values (i.e. as close to plus or minus one as possible). Numerical algorithms will fail if you feed them crazy values. – Erwin Kalvelagen Oct 09 '17 at 01:14
  • Many thanks. So I divide the values in the original objective function by 10e13 and get c = [-0.7890424934354171875,-0.7890424934354274414,-0.7890424934‌​354246093], still I couldn't get a correct answer. I think the reason is that the precision is lost during the computation of ILP. Might you please give me some advice if I have to deal with the case that coefficients are crazy or very close to each other like in the example where high precision is required? @Erwin Kalvelagen – Sue Oct 09 '17 at 02:00
  • 1
    Every general solver will have trouble with your scaled c as internally doubles or maybe some tiny extensions (80 bits? 128? at least for solving some equations) are used. Apart from looking for specialized solvers (rational arithmetic might be some approach; not really recommended), the problem of yours is more linked to modelling than solving. – sascha Oct 09 '17 at 02:05
  • 1
    Replace `min -7890424934354.171875*x1 -7890424934354.274414*x2 -7890424934354.246093*x3` by `min -7890424934354 -0.171875*x1 -0.274414*x2 -0.246093*x3` and drop the constant term. (We can do this because `x1 + x2 + x3 = 1`). – Erwin Kalvelagen Oct 09 '17 at 02:08
  • Yes, this is a good idea. I have also truncated the values the same way since x1 + x2 + x3 = 1 in the example. I was wondering if there is a general way or some other tools to solve this kind of problem since the constraint may be something else besides x1 + x2 + x3 = 1. @Erwin Kalvelagen – Sue Oct 09 '17 at 05:43
  • No, these original objective coefficients are essentially the same. The solver cannot distinguish them. Remember that besides floating point representation and arithmetic issues, numeric algorithms deal with tolerances. For normal practical models, this is usually not an issue. This looks like a somewhat artificial problem. – Erwin Kalvelagen Oct 09 '17 at 05:53
  • Thank you for your advice. I should think more about my original model instead of trying to find specialized solvers which support high precision. @sascha – Sue Oct 09 '17 at 05:59
  • Oh, got it. I happened to get the values in the objective function, which is not the normal case. Thank you very much. @Erwin Kalvelagen – Sue Oct 09 '17 at 06:05

0 Answers0