[Help-glpk] Re: Simplex method infeasible, IP method feasible

[Andrew Makhorin]

A simple analysis shows that your problem being primal infeasible is
nevertheless almost feasible.

Glpk lp presolver does not recover non-optimal and infeasible solutions,
so disable it (in glpsol it is enabled by default):

   glpsol --lpt infeas.lpt --nopresol -o infeas.sol

and look at the solution reported by the simplex solver:

========================================================================
Problem:    PROBLEM
Rows:       6
Columns:    3
Non-zeros:  12
Status:     INFEASIBLE (FINAL)
Objective:  obj = 68.85504317 (MINimum)

  Row name   St   Activity     Lower bound   Upper bound    Marginal
------------ -- ------------- ------------- ------------- -------------
r(1)         NU             0                           0            -1
r(2)         B         68.855                      68.855 
r(3)         B        -43.537                     39.4829 
r(4)         B       -93.4892                     5.23378 
r(5)         NU      -73.0421                    -73.0421      -1.44146
r(6)         NU       76.0438                     76.0438     -0.479097

Column name  St   Activity     Lower bound   Upper bound    Marginal
------------ -- ------------- ------------- ------------- -------------
x(1)         B         68.855        -1e+06         1e+06 
x(2)         B        92.4179        -1e+06         1e+06 
x(3)         B       -23.5629        -1e+06         1e+06 

Karush-Kuhn-Tucker optimality conditions:

KKT.PE: max.abs.err. = 1.42e-14 on row 2
        max.rel.err. = 4.49e-15 on row 1
        High quality

KKT.PB: max.abs.err. = 8.46e-06 on row 2
        max.rel.err. = 1.21e-07 on row 2
        Medium quality

KKT.DE: max.abs.err. = 1.81e-16 on column 3
        max.rel.err. = 1.81e-16 on column 3
        High quality

KKT.DB: max.abs.err. = 0.00e+00 on row 0
        max.rel.err. = 0.00e+00 on row 0
        High quality

End of output
========================================================================

As it seen from KKT optimality conditions some bound constraints cannot
be satisfied within a given tolerance, and row 2 has the largest (among
other rows and columns) relative error 1.21e-7. (The difference between
primal value of r(2) and its upper bound which is violated is not seen
because only 6 decimal places are printed.) If all residuals were less
than tol_bnd (1e-7 by default), the simplex solver would say that the
solution is optimal. The floating-point arithmetic is inexact, so if you
agree that 1.21e-7 is the same as zero, your problem is primal feasible
and therefore the solution above is optimal (because other optimality
conditions are not violated). If you do not agree, your problem has no
feasible solutions. To agree or not to agree depends on the nature of
your problem, in particular, on the accuracy of data.

Unfortunately glpk has no feature to perform infeasibility analysis
(I hope to implement that in the future), but it is easy to do that
manually. Any infeasible problem can be made feasible by introducing
non-negative artificial variables, which play the role of residuals,
by one for each row (constraint). In any primal feasible solution of the
original problem these residuals must be zero, so we can understand what
causes infeasibility minimizing the sum of such artificial variables:

========================================================================
\* Problem: Unknown *\

Minimize
 obj: y1 + y2 + y3 + y4 + y5 + y6

Subject To
 r(1): -y1 - x(1) + x(2) + x(3) <= 0
 r(2): -y2 + x(1) <= 68.855034718
 r(3): -y3 - 0.222520933956 x(2) + 0.974927912182 x(3) <= 39.4828617494
 r(4): -y4 - 0.900968867902 x(2) + 0.433883739118 x(3) <= 5.23377925511
 r(5): -y5 - 0.900968867902 x(2) - 0.433883739118 x(3) <= -73.0421245614
 r(6): -y6 + 0.623489801859 x(2) - 0.781831482468 x(3) <= 76.0438443135

Bounds
 -1000000 <= x(1) <= 1000000
 -1000000 <= x(2) <= 1000000
 -1000000 <= x(3) <= 1000000

End
========================================================================

Solving the auxiliary problem we have:

lpt_read_prob: 9 variables, 6 constraints
lpt_read_prob: 19 lines were read
gm_scal: max / min = 4.494e+00
gm_scal: max / min = 3.016e+00
lpx_adv_basis: size of triangular part = 6
      0:   objval =  -7.922441352e+05   infeas =   1.000000000e+00 (0)
      5:   objval =   4.008567879e+01   infeas =   0.000000000e+00 (0)
*     5:   objval =   4.008567879e+01   infeas =   0.000000000e+00 (0)
*     7:   objval =   5.865824880e-06   infeas =   0.000000000e+00 (0)
OPTIMAL SOLUTION FOUND

The minimal sum of residuals is non-zero (i.e. some y[i] are non-zero),
ergo the original problem has no primal feasible solution.

Now look at the optimal solution of this auxiliary problem:

========================================================================
Problem:    PROBLEM
Rows:       6
Columns:    9
Non-zeros:  18
Status:     OPTIMAL
Objective:  obj = 5.86582488e-06 (MINimum)

  Row name   St   Activity     Lower bound   Upper bound    Marginal
------------ -- ------------- ------------- ------------- -------------
r(1)         NU             0                           0      -0.69374
r(2)         NU        68.855                      68.855      -0.69374
r(3)         B        -43.537                     39.4829 
r(4)         B       -93.4892                     5.23378 
r(5)         NU      -73.0421                    -73.0421            -1
r(6)         NU       76.0438                     76.0438     -0.332369

Column name  St   Activity     Lower bound   Upper bound    Marginal
------------ -- ------------- ------------- ------------- -------------
y1           NL             0             0                     0.30626
y2           NL             0             0                     0.30626
y3           NL             0             0                           1
y4           NL             0             0                           1
y5           B    5.86582e-06             0               
y6           NL             0             0                    0.667631
x(1)         B         68.855        -1e+06         1e+06 
x(2)         B        92.4179        -1e+06         1e+06 
x(3)         B       -23.5629        -1e+06         1e+06 
========================================================================

We see that four rows, namely, r(1), r(2), r(5), r(6), are active, and
they all have negative reduced costs of right-hand sides. This means
that it is sufficient to slightly increase the right-hand side of any of
these four rows, and the original problem will become primal feasible.
For example, if we increase the rhs of r(2) by 1, the sum of residuals
is decreased by 0.69374. Since we want to make the sum being zero, it is
sufficient, for example, to increase the right-hand side of r(2) at
least by 5.86e-6 / 0.69374 = 8.45e-6. The rhs of r(2) is 68.855034718,
and 8.45e-6 is relatively very small quantity which allows eliminating
primal infeasibility. Try to replace the rhs of r(2) by 68.856, and you
will see that the simplex will report you optimal solution. That is in
what sense your problem is almost primal feasible.

The glpk interior point solver (as well as cplex) reports the optimal
solution due to happy chance, because as a rule ip methods are more
tolerant to primal infeasibility. In principle, it is sufficient to
change the corresponding tolerance or make the constraints in your
problem a bit tighter, and the ip solver also will tell you that your
problem is infeasible.

Hope this helps you to understand what is happening :+)

Best regards,

Andrew Makhorin


----- Original Message ----- 
From: Nicolo' Giorgetti <address@hidden>
To: <address@hidden>
Sent: Friday, September 12, 2003 7:35 PM
Subject: [Help-glpk] Re: Simplex method infeasible, IP method feasible


> On Fri, 12 Sep 2003 14:41:32 +0400 Andrew Makhorin <address@hidden> wrote:
> 
> > From: Nicolo' Giorgetti <address@hidden>
> > To: <address@hidden>
> > Subject: Re: [Help-glpk] Re: Simplex method infeasible, IP method
> > feasible Date: Friday, September 12, 2003 2:20 PM
> > 
> > >I'm doing some errors in defining the problem but I'm not able
> > >to see where they are. Do you have any suggestion ?
> > >Thank you very much.
> > 
> > CONFER!
> 
> Dear Andrew,
> 
> Sorry, I have another problem similar to the previous one.
> 
> Could you solve it both with the simplex method and the interior point 
> method ? 
> 
> --------- problem: output_lpt.lpt -------------------
> 
> \* Problem: Unknown *\
> 
> Minimize
>  obj: + x(1)
> 
> Subject To
>  r(1): - x(1) + x(2) + x(3) <= 0
>  r(2): + x(1) <= 68.855034718
>  r(3): - 0.222520933956 x(2) + 0.974927912182 x(3) <= 39.4828617494
>  r(4): - 0.900968867902 x(2) + 0.433883739118 x(3) <= 5.23377925511
>  r(5): - 0.900968867902 x(2) - 0.433883739118 x(3) <= -73.0421245614
>  r(6): + 0.623489801859 x(2) - 0.781831482468 x(3) <= 76.0438443135
> 
> Bounds
>  -1000000 <= x(1) <= 1000000
>  -1000000 <= x(2) <= 1000000
>  -1000000 <= x(3) <= 1000000
> 
> End
> 
> -------------------------------------------------------
> 
> 
> I solved it using glpsol and I obtained two different outputs:
> 
> --------------- Simplex ------------------------
> $glpsol --lpt output_lpt.lpt
> 
> lpt_read_prob: reading LP data from `output_lpt.lpt'...
> lpt_read_prob: 3 variables, 6 constraints
> lpt_read_prob: 19 lines were read
> lpx_simplex: original LP has 6 rows, 3 columns, 12 non-zeros
> lpx_simplex: presolved LP has 5 rows, 3 columns, 11 non-zeros
> gm_scal: max / min = 4.494e+000
> gm_scal: max / min = 3.016e+000
> lpx_adv_basis: size of triangular part = 5
>       0:   objval = -2.000000000e+006   infeas =  1.000000000e+000 (0)
>       4:   objval =  6.885503472e+001   infeas =  2.000844258e-012 (0)
> *     4:   objval =  6.885503472e+001   infeas =  5.401980374e-006 (0)
> *     5:   objval =  6.885504317e+001   infeas =  7.015617570e-006 (0)
> lpx_simplex: numerical instability (primal simplex, phase II)
>       5:   objval =  6.885504317e+001   infeas =  1.000000000e+000 (0)
>       5:   objval =  6.885504317e+001   infeas =  1.000000000e+000 (0)
> PROBLEM HAS NO FEASIBLE SOLUTION
> lpx_simplex: cannot recover undefined or non-optimal solution
> Time used:   0.0 secs
> Memory used: 0.1M (80744 bytes)
> 
> ------- Interior point method ------------------------
> 
> $ glpsol --interior --lpt output_lpt.lpt
> 
> lpt_read_prob: reading LP data from `output_lpt.lpt'...
> lpt_read_prob: 3 variables, 6 constraints
> lpt_read_prob: 19 lines were read
> lpx_interior: original LP problem has 6 rows and 3 columns
> lpx_interior: transformed LP problem has 9 rows and 12 columns
> Factorizing S = A*A' symbolically...
> nz(A) = 24
> nz(S) = 32 (upper triangular part)
> nz(U) = 32
> Guessing initial point...
> Optimization begins...
>   0: F =  2.522900440e+005; rpi = 3.0e-001; rdi = 1.1e+000; gap =
> 2.0e-001  
>   1: F = -4.073474860e+004; rpi = 9.2e-002; rdi = 1.1e-001; gap
> = 1.8e-001  
>   2: F = -1.149990823e+003; rpi = 1.1e-002; rdi = 1.3e-002;
> gap = 2.4e-002  
>   3: F = -4.535524941e+001; rpi = 1.1e-003; rdi = 1.3e-003; gap =
> 2.5e-003   
>   4: F =  6.094629906e+001; rpi = 1.2e-004; rdi= 1.3e-004; gap =
> 3.2e-004  
>   5: F =  7.301168842e+001; rpi = 1.5e-005; rdi = 1.3e-005; gap =
> 5.8e-005  
>   6: F =  7.071401526e+001; rpi = 2.3e-006; rdi = 1.7e-006; gap =
> 7.9e-006  
>   7: F =  6.904379861e+001; rpi= 2.3e-007; rdi = 1.7e-007; gap= 8.1e-007
>   8: F =  6.887392223e+001; rpi = 2.3e-008; rdi= 1.7e-008; gap= 8.1e-008
>   9: F = 6.885693739e+001; rpi = 2.3e-009; rdi= 1.7e-009; gap= 8.1e-009
> OPTIMAL SOLUTION FOUND 
> Time used:   0.0 secs 
> Memory used: 0.1M(58588 bytes)
> 
> 
> I solved it with the barrier algorithm of CPLEX and I have obtained the
> following output:
> 
> ------------ from CPLEX --------------------
> 
> Problem 'output_lpt.lpt' read.
> Read time =    0.10 sec.
> Tried aggregator 1 time.
> LP Presolve eliminated 1 rows and 0 columns.
> Reduced LP has 5 rows, 3 columns, and 11 nonzeros.
> Presolve time =    0.02 sec.
> Number of nonzeros in lower triangle of A*A' = 10
> Using Approximate Minimum Degree ordering
> Total time for automatic ordering = 0.00 sec.
> Summary statistics for Cholesky factor:
>   Rows in Factor            = 5
>   Integer space required    = 5
>   Total non-zeros in factor = 15
>   Total FP ops to factor    = 55
>  Itn      Primal Obj        Dual Obj  Prim Inf Upper Inf  Dual Inf      
>       0 -3.9574734e+005 -1.1000138e+007 3.10e+006 2.08e+005 1.00e+001
>    1 -3.1615429e+005 -1.7627440e+006 7.28e+005 4.89e+004 1.98e+000
>    2  2.3392532e+003 -9.5588211e+004 3.41e+004 2.29e+003 6.48e-015
>    3  6.9185003e+001 -1.8265269e+002 5.61e+000 3.77e-001 6.23e-015
>    4  6.8912782e+001  6.8129432e+001 1.09e+000 7.31e-002 4.27e-015
>    5  6.8846880e+001  6.8615031e+001 3.32e-002 2.13e-003 1.67e-013
>    6  6.8854183e+001  6.8835404e+001 1.96e-003 2.26e-015 4.06e-015
>    7  6.8854957e+001  6.8853283e+001 1.86e-004 2.33e-010 1.96e-015
>    8  6.8855028e+001  6.8854909e+001 2.43e-005 1.63e-016 1.24e-015
>    9  6.8855031e+001  6.8855056e+001 1.89e-005 3.09e-015 2.82e-015
>   10  6.8855035e+001  6.8855067e+001 9.21e-006 3.58e-015 1.90e-015
>   11  6.8855035e+001  6.8855071e+001 7.90e-006 6.67e-015 1.77e-015
>  Barrier cannot determine infeasibility.
> Barrier time =    0.05 sec.
> 
> Primal crossover.
>   Primal:  Fixed no variables.
>   Dual:  Fixing 1 variable.
>         0 DMoves:  Infeasibility 0.00000000e+000  Objective
> 6.88550432e+001  Dual:  Pushed 1, exchanged 0.
> Using devex.
> Removing shift (1).
> Total crossover time =    0.08 sec.
> 
> Primal simplex - Infeasible:  Infeasibility =   8.4553621207e-006
> Solution time =    0.18 sec.  Iterations = 0 (0)
> 
> 
> Best Regards,
> Nicolo'.
> 
> 
> _______________________________________________
> Help-glpk mailing list
> address@hidden
> http://mail.gnu.org/mailman/listinfo/help-glpk
>