Re: Solver delivers wrong answer when 2 constraints are close

Hello João !

If the model/data that you tested could be public, do you mind share then ?

Cheers !

On 4/3/21 19:00, João Flávio de Freitas Almeida wrote:
> Dear all,
>
> Maybe it is a design flaw, or an AMPL issue but recently, I run the 
> same model with Glpk and Ampl-Gurobi.
> The optimal solution were different.
>
> *GLPK: *3.328774238e+09
> *AMPL-Gurobi: *3.360308548e+09
>
> # ========================================================== #
> GLPSOL: GLPK LP/MIP Solver, v4.65
> Parameter(s) specified in the command line:
>  -m v9_3flass_glpk.mod -d v9_dados.dat -d 12_dados.dat -d 
> 12_distancias.txt
>  --cuts --tmlim 7200 --log v9_log.log
> Reading model section from v9_3flass_glpk.mod...
> 674 lines were read
> Reading data section from v9_dados.dat...
> 112 lines were read
> Reading data section from 12_dados.dat...
> 61 lines were read
> Reading data section from 12_distancias.txt...
> 491 lines were read
> Checking (line 184)...
> Checking (line 242)...
> Checking (line 243)...
> Checking (line 244)...
> Generating of...
> Generating R1a...
> Generating R2...
> Generating R3...
> Generating R4...
> Generating R5...
> Generating R6...
> Generating R7a...
> Generating R7b...
> Generating R7c...
> Generating R8...
> Generating R9...
> Generating R10...
> Generating R10a...
> Generating R11...
> Generating R11a...
> Generating R11b...
> Generating R12...
> Generating R12b...
> Generating R13...
> Generating R14...
> Generating R15...
> Generating R16...
> Generating R17...
> Generating R18...
> Generating R19...
> Generating R20...
> Generating R21...
> Model has been successfully generated
> GLPK Integer Optimizer, v4.65
> 5139 rows, 1248 columns, 23444 non-zeros
> 608 integer variables, all of which are binary
> Preprocessing...
> 165 constraint coefficient(s) were reduced
> 766 rows, 1019 columns, 15790 non-zeros
> 509 integer variables, all of which are binary
> Scaling...
>  A: min|aij| =  1.000e-02  max|aij| =  2.006e+01  ratio =  2.006e+03
> GM: min|aij| =  7.199e-01  max|aij| =  1.389e+00  ratio =  1.930e+00
> EQ: min|aij| =  5.209e-01  max|aij| =  1.000e+00  ratio =  1.920e+00
> 2N: min|aij| =  2.500e-01  max|aij| =  1.680e+00  ratio =  6.720e+00
> Constructing initial basis...
> Size of triangular part is 766
> Solving LP relaxation...
> GLPK Simplex Optimizer, v4.65
> 766 rows, 1019 columns, 15790 non-zeros
>       0: obj =   2.904848487e+09 inf =   3.040e+02 (75)
>      79: obj =   3.630095622e+09 inf =   4.441e-16 (0)
> *   130: obj =   3.328774238e+09 inf =   7.881e-15 (0)
> OPTIMAL LP SOLUTION FOUND
> Integer optimization begins...
> Long-step dual simplex will be used
> Gomory's cuts enabled
> MIR cuts enabled
> Cover cuts enabled
> Number of 0-1 knapsack inequalities = 95
> Clique cuts enabled
> Constructing conflict graph...
> Conflict graph has 182 + 70 = 252 vertices
> +   130: mip =     not found yet >=              -inf    (1; 0)
> Cuts on level 0: gmi = 2; mir = 2;
> Cuts on level 10: gmi = 2; mir = 2;
> +   161: >>>>>   3.328774238e+09 >= 3.328774238e+09   0.0% (11; 0)
> +   161: mip =   3.328774238e+09 >=     tree is empty 0.0% (0; 21)
> INTEGER OPTIMAL SOLUTION FOUND
> Time used:   0.3 secs
> Memory used: 12.2 Mb (12789041 bytes)
> Model has been successfully processed
>
> # ========================================================== #
>
> *Gurobi 9.1.1 (win64) logging started Thu Mar  4 14:41:08 2021*
>
> Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)
> Thread count: 2 physical cores, 4 logical processors, using up to 4 
> threads
> Optimize a model with 766 rows, 1021 columns and 15790 nonzeros
> Variable types: 510 continuous, 511 integer (511 binary)
> Coefficient statistics:
>   Matrix range     [1e-02, 2e+01]
>   Objective range  [6e+05, 4e+08]
>   Bounds range     [1e+00, 1e+00]
>   RHS range        [4e-02, 2e+01]
> Found heuristic solution: objective 3.666169e+09
> Presolve removed 759 rows and 1012 columns
> Presolve time: 0.03s
> Presolved: 7 rows, 9 columns, 18 nonzeros
> Found heuristic solution: objective 3.360309e+09
> Variable types: 0 continuous, 9 integer (9 binary)
>
> Root relaxation: cutoff, 2 iterations, 0.01 seconds
>
>     Nodes    |    Current Node    |     Objective Bounds  |     Work
>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd Gap | It/Node Time
>
>      0     0     cutoff    0      3.3603e+09 3.3603e+09  0.00%     -    0s
>
> Explored 0 nodes (2 simplex iterations) in 0.09 seconds
> Thread count was 4 (of 4 available processors)
>
> Solution count 2: 3.36031e+09 3.66617e+09
>
> Optimal solution found (tolerance 1.00e-04)
> Best objective 3.360308547685e+09, best bound 3.360308547685e+09, gap 
> 0.0000%
> Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)
> Thread count: 2 physical cores, 4 logical processors, using up to 4 
> threads
> Optimize a model with 766 rows, 1021 columns and 15790 nonzeros
> Coefficient statistics:
>   Matrix range     [1e-02, 2e+01]
>   Objective range  [6e+05, 4e+08]
>   Bounds range     [1e+00, 1e+00]
>   RHS range        [4e-02, 2e+01]
> Iteration    Objective       Primal Inf.    Dual Inf.  Time
>        0    3.3603085e+09   2.034244e+01   0.000000e+00  0s
>       45    3.3603085e+09   0.000000e+00   0.000000e+00  0s
>
> Solved in 45 iterations and 0.01 seconds
> *
> *
> Optimal objective  3.360308548e+09
>
> ========================================================
>
> Bests,
> --
> João Flávio de Freitas Almeida
>
>
> Em qui., 4 de mar. de 2021 às 14:20, Thiago Neves 
> <[email protected] <mailto:[email protected]>> escreveu:
>
>     Thank you so much for the explanation Antti.
>     You are right, this seems more like a "design flaw" than a "bug".
>     I'm sorry for the misconception.
>
>     Why does "eps" need to be so big (1e-3)? Is there a way to set a
>     fixed threshold in the command line (glpsol) or at least the 1e-3
>     part?
>
>     In the case (1), wouldn't it be better to choose max(l[q], l'[q])?
>
>     Att,
>
>     *Thiago H. Neves*
>     (31) 98608-0666
>
>
>
>     ---------- Forwarded message ---------
>     De: *Antti Lehtila* <[email protected]
>     <mailto:[email protected]>>
>     Date: qui., 4 de mar. de 2021 às 12:44
>     Subject: Re: Solver delivers wrong answer when 2 constraints are close
>     To: <[email protected] <mailto:[email protected]>>
>
>
>     Hi,
>
>     I think it works fully as documented, and so *per design*. Singleton
>     rows are treated as column bounds by the preprocessor. See
>     documentation
>     for *npp_implied_lower*:
>     *---------------
>     *  Processing implied column lower bound l'[q] includes the following
>     *  cases:
>     *
>     *  1) if l'[q] < l[q] + eps, implied lower bound is redundant;
>     *
>     *  2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound
>     *     l[q] can be strengthened by replacing it with l'[q]. If in this
>     *     case new column lower bound becomes close to current column
>     upper
>     *     bound u[q], the column can be fixed on its upper bound;
>     *
>     *  3) if l'[q] > u[q] + eps, implied lower bound violates current
>     *     column upper bound u[q], in which case the problem has no primal
>     *     feasible solution. */
>     *---------------
>     The lower bound can have only a single value, but if you define
>     multiple
>     values for a column lower bound, they must of course be processed in
>     some order.  In this case, the lower bound is first defined
>     l(q)=0, then
>     l'(q)=1, and finally l'(q)=1.0001. The third value for the bound
>     is thus
>     considered *redundant*, as per design, and so the second value l(q)=1
>     remains in effect. This is because *eps* is defined as 1e-3 + 1e-6 *
>     fabs(q->lb)).
>
>     I guess it may be a considered a design flaw, but I think it
>     should not
>     be called a bug, as it is working as designed and documented. 
>     Besides,
>     I think one should use the Bounds section for bounds, instead of
>     using
>     multiple constraints for defining a single lower bound.
>
>     Best,
>     Antti
>     _______________________
>     On 04.03.2021 11:38, Domingo Alvarez Duarte wrote:
>     >
>     > Testing this problem I discover that if we change the order of
>     > constraint declarations it seems to give the expected answer as
>     stated
>     > by Thiago (what I think could be another bug).
>     >
>     > ====
>     >
>     > /param min_bound default 0;/
>     > /var x >= 0;
>     >
>     > minimize y: x;/
>     > *
>     > */*s.t. PART_MIN_X: x >= 1 + min_bound;*
>     > /
>     > /*s.t. LIM_INF_X: x >= 1;
>     > *
>     > /
>     > /solve;
>     > display min_bound;
>     > display x; # EXPECTED RESULT: X ==  1.0001
>     >
>     > data;
>     > param min_bound := 1e-4;
>     > end;/
>     >
>     > ====
>     >
>     > Output:
>     >
>     > ====
>     >
>     > x.val = 1.0001
>     >
>     > ====
>     >
>     > On 3/3/21 19:19, Thiago Neves wrote:
>     >> Hi.
>     >> I've found a strange behaviour in glpk which I don't know how
>     to fix
>     >> nor how to contour it. It seems like GLPK can't distinguish
>     >> constraints that differs from about 1e-4.
>     >>
>     >> Follows simple examples that explain and reproduce the problem.**
>     >> *
>     >> *
>     >> *The first model gives the desired answer (x = 1.0001):*
>     >> /
>     >> param min_bound default 0;/
>     >> /var x >= 0;
>     >>
>     >> minimize y: x;/
>     >> /*
>     >> s.t. PART_MIN_X: x >= 1 + min_bound;*
>     >>
>     >> solve;
>     >> display min_bound;
>     >> display x; # EXPECTED RESULT: X ==  1.0001
>     >>
>     >> data;
>     >> param min_bound := 1e-4;
>     >> end;
>     >> /
>     >> /_____________________________________/
>     >> /OUTPUT:/
>     >> /x.val = 1.0001/
>     >> /_____________________________________ /
>     >>
>     >> *Now, if I add a second constraint "close" to the first one, the
>     >> solver will deliver an answer that is actually infeasible:*
>     >>
>     >> /param min_bound default 0;/
>     >> /var x >= 0;
>     >>
>     >> minimize y: x;/
>     >>
>     >> *s.t. LIM_INF_X: x >= 1;
>     >>
>     >> */*s.t. PART_MIN_X: x >= 1 + min_bound;*
>     >>
>     >> solve;
>     >> display min_bound;
>     >> display x; # EXPECTED RESULT: X ==  1.0001
>     >>
>     >> data;
>     >> param min_bound := 1e-4;
>     >> end;/
>     >> /_____________________________________/
>     >> /OUTPUT:/
>     >> x.val = 1
>     >> /_____________________________________ /
>     >>
>     >> *If I change the "min_bound" parameter to 1e-2, the second model
>     >> works as expected (x = 1.01):*
>     >>
>     >> /param min_bound default 0;/
>     >> /
>     >> /
>     >> /var x >= 0;
>     >>
>     >> minimize y: x;/
>     >> *
>     >> *
>     >> *s.t. LIM_INF_X: x >= 1;
>     >>
>     >> */*s.t. PART_MIN_X: x >= 1 + min_bound;*
>     >>
>     >> solve;/
>     >> /
>     >> display x; # EXPECTED RESULT: X ==  1.01
>     >>
>     >> data;
>     >> param min_bound := 1e-2;
>     >> end;/
>     >> /_____________________________________/
>     >> /OUTPUT:/
>     >> x.val = 1.01
>     >> /_____________________________________ /
>     >> /
>     >> /
>     >>
>     >> Att,
>     >>
>     >> *Thiago H. Neves*
>     >> (31) 98608-0666
>     >>
>
>