Re: Solver delivers wrong answer when 2 constraints are close

[Andrew Makhorin]

On Thu, 2021-03-04 at 17:47 +0100, Heinrich Schuchardt wrote:
> On 3/4/21 4:44 PM, Antti Lehtila wrote:
> > 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
> 
> In this line an apostrophe is missing.
> 
>    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;
> 
> It is this strengthening that fails:
> 
> src/npp/npp3.c:567
> eps = (q->is_int ? 1e-3 : 1e-3 + 1e-8 * fabs(q->lb));
> 
> Set eps to 1E-5 and you are fine.
> 
> Or run with --nopresol.
> 
> @Andrew:
> Shouldn't 1E-5 be good enough?
> Why do we need eps > 0?


Both 1e-3 and 1e-5 are good enough.

If the result is sensitive to such small deviations, this only tells
that the model is badly formulated.


> 
> Best regards.
> 
> Heinrich
> 
> > *
> > *  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
> > > > 
> > 
> > 
> 
>