scipopt / papilo

Parallel Presolve for Integer and Linear Optimization
GNU Lesser General Public License v3.0
64 stars 17 forks source link

The permutation of constraints and variables can cause infeasible. #29

Closed liangbug closed 1 year ago

liangbug commented 1 year ago

I am dealing with a problem on SCIP 8.0.3, but it can't get the optimal value, which turns out infeasible. I permutate constraints and variables, then get the optimal value. What may cause it to happen? Are there any suggestions to prevent this error?

ambros-gleixner commented 1 year ago

This should of course not happen. Please file a bug report with the exact version and settings you are using, either here or via https://scipopt.org/bugs.php.

liangbug commented 1 year ago

Unfortunately the data is confidential, I’ll try to provide some simulated data. The infeasibility occurs in Papilo propagation presolver, SingleRow.hpp:102, the method checkStatus() first kInfeasible.

version: SCIP 8.0.3, PaPILO 2.1.2

setting: presolving/donotmulagger=True

log: fail.log success.log

alexhoen commented 1 year ago

From the log it is nearly impossible to judge what causes the problem. This is the time this appears, so providing simulated data (or a subset of the problem) would be great to fix this problem.

alexhoen commented 1 year ago

Additional: The SingleRow identifies problems as infeasible if it contains only a variable and the row forces a value that is not within its domain defined by the bounds. An option that we can also offer you is to explain to you how to debug the code and then we can hopefully fix the problem with the information gathered. In this case, it is necessary to track all the reductions performed on this constraint and identify the culprit reduction. (This is the standard way I would start and usually it works, but can not make any promise...) If that is an option for you I will send you advice on how to do this. It should not be that hard.

liangbug commented 1 year ago

It sounds ok to understand how to debug the code with your help. Please kindly show me how to proceed. Thank you.

alexhoen commented 1 year ago

So I think there are three steps, which I will explain first briefly and then go into detail:

  1. Get the (original) index of the infeasible single row.
  2. Enable the verbose PaPILO output to print all the reductions performed on the instance.
  3. Search for the row in the output and check which reduction is falsely applied.

To make things a little more easy would recommend the following two things:

  1. Check if this bug also appears if only PaPILO is used. (https://github.com/scipopt/papilo) Since PaPILO can not read cip files it might be necessary to convert the instance to a mps file. Use the command-line SCIP to read the instance and then write it out as mps file. (build/bin/papilo presolve -f INSTANCE -p settings/default.set)
  2. Disabling presolvers that do not contribute to the bug reduces the log and the potential false reductions. To do so disable as many presolvers as possible and check if the bug still appears. For the PaPILO settings this is PRESOLVERNAME.enabled = 0 (colsingleton,coefftightening, propagation, simpleprobing, parallelrows, stuffing, dualfix, fixcontinuous, simplifyineq, doubletoneq, dualinfer, implint, domcol, probing, substitution, sparsify). For SCIP the settings would be: presolving/milp/x with x = {enableprobing, enabledomcol, enableparallelrows, enabledualinfer, enablemultiaggr}

So now I would explain the three necessary steps:

  1. Set a breakpoint in SingleRow:102. Move on debug-frame up and there you can should see and integer variable r or row. PaPILO shrinks the matrix so this variable contains the row index only in the reduced matrix. To obtain the original row value please evaluate problem.postsolve.origrow_mapping[row]. (maybe it is worth taking a closer look at the original row: problem.postsolve.getOriginalProblem().getConstraintMatrix().getRowCoefficients()[originalrow])
  2. now run PaPILO with verbose settings by setting the output to 4. This results in a long list of numbers. To convert it something human readable please use the parser in PaPILO: https://github.com/scipopt/papilo/tree/master/check/parser . Save the log output to test.out and then replace the main in LogFileParser with
    if __name__ == '__main__':
    (applied_tsx, dependency_matrix, canceled_tsx, inactive_tsx, errors, rounds, rounds_per_presolver,
    same_appeareances_per_presolver, conflict_appearances_per_presolver) = run_analysis_for_file("test.out", True)

    This should result in something that is hopefully understandable.

  3. Now search for the appearances of the original_row in the log file. And maybe it would be wise to search also for the variable indexes in the log file as well. if you reached this point we should definitely discuss the results.

I hope you can follow the steps. if there are any problem do not hestitate the ask.

liangbug commented 1 year ago

In my original problem, I found that

  1. the continuous variable cannot be fixed to a constant if the difference between the upper and lower bounds is less than epsilon after several rounds of presolving.

  2. In the Papilo Simple Probing Presolver. If the last variable of the constraint is continuous and the difference between upper and lower bounds of the last variable is less than epsilon, the constraint is reduced to a singleton row. Then the binary variable of the singleton row is fixed to 0, since the rhs of the singleton row is less than epsilon. But the constraint shouldn't be reduced to a singleton row.

I can't reproduce the [1.] problem easily. And I created the following simple model to show the [2.] error.

# version: SCIP 8.0.3, Papilo 2.1.2
from pyscipopt import Model, SCIP_PARAMSETTING

if __name__ == '__main__':
    model = Model()

    model.setIntParam('display/verblevel', 5)
    model.setHeuristics(SCIP_PARAMSETTING.OFF)
    model.setPresolve(SCIP_PARAMSETTING.OFF)
    model.setIntParam('presolving/milp/maxrounds',1)
    model.setIntParam('presolving/maxrounds', 1)

    x = model.addVar(name='x', vtype='B')
    y = model.addVar(name='y', vtype='C', lb=0-1e-12, ub=2+1e-12)
    z = model.addVar(name='z', vtype='C', lb=1-1e-12, ub=1+1e-12)

    model.addCons(2*x + y + z == 3, name='c1')
    model.addCons(3*x + 2*y + z >= 4 , name='c2')

    model.setObjective(x, 'maximize')

    model.optimize()
    model.printBestSol()
    # x should be 1, but 0.
    # x would be fixed to 0 after Simple Probling.

log

LP Solver <Soplex 6.0.3>: barrier convergence tolerance cannot be set -- tolerance of SCIP and LP solver may differ
LP Solver <Soplex 6.0.3>: fastmip setting not available -- SCIP parameter has no effect
LP Solver <Soplex 6.0.3>: number of threads settings not available -- SCIP parameter has no effect
transformed problem has 3 variables (1 bin, 0 int, 0 impl, 2 cont) and 2 constraints
      2 constraints of type <linear>

original problem has 6 active (100%) nonzeros and 6 (100%) check nonzeros

presolving:
   (0.0s) running MILP presolver

starting presolve of problem model:
  rows:     2
  columns:  3
  int. columns:  1
  cont. columns:  2
  nonzeros: 6

round 0   ( Trivial  ):    0 del cols,    0 del rows,    0 chg bounds,    0 chg sides,    0 chg coeffs,    0 tsx applied,    0 tsx conflicts
Presolver simpleprobing applying 
row -10 col 1 val -2.0
row -1 col 0 val 2.0
tsx
modified rows: 0,1,
row -10 col 2 val -1.999511667349907e-13
row -1 col 0 val 1.0000000000001
tsx
modified rows: 0,1,
round 1   (  Medium  ):    3 del cols,    2 del rows,    1 chg bounds,    8 chg sides,    8 chg coeffs,    2 tsx applied,    0 tsx conflicts
presolved 2 rounds:    3 del cols,    2 del rows,    1 chg bounds,    8 chg sides,    8 chg coeffs,    2 tsx applied,    0 tsx conflicts

          presolver     nb calls   success calls(%)    nb transactions     tsx applied(%)  execution time(s) 
       colsingleton            1                0.0                  0                0.0              0.000
    coefftightening            1                0.0                  0                0.0              0.000
        propagation            1                0.0                  0                0.0              0.000
      simpleprobing            1              100.0                  2              100.0              0.000
       parallelrows            1                0.0                  0                0.0              0.000
           stuffing            1                0.0                  0                0.0              0.000
            dualfix            1                0.0                  0                0.0              0.000
      fixcontinuous            0                0.0                  0                0.0              0.000
       simplifyineq            0                0.0                  0                0.0              0.000
        doubletoneq            1                0.0                  0                0.0              0.000
            implint            0                0.0                  0                0.0              0.000
          dualinfer            0                0.0                  0                0.0              0.000
            probing            0                0.0                  0                0.0              0.000
             domcol            0                0.0                  0                0.0              0.000
       substitution            0                0.0                  0                0.0              0.000

reduced problem:
  reduced rows:     0
  reduced columns:  0
  reduced int. columns:  0
  reduced cont. columns:  0
  reduced nonzeros: 0
Solution passed validation
problem is solved [optimal solution found] [objective value: 0.0 (double precision)]
   (0.0s) MILP presolver (2 rounds): 2 aggregations, 1 fixings, 0 bound changes

presolved problem has 0 active (0%) nonzeros and 0 (0%) check nonzeros

presolving (1 rounds: 1 fast, 1 medium, 1 exhaustive):
 3 deleted vars, 2 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 0 cliques
transformed 1/1 original solutions to the transformed problem space
Presolving Time: 0.00

SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes      : 0
Primal Bound       : +0.00000000000000e+00 (1 solutions)
Dual Bound         : +0.00000000000000e+00
Gap                : 0.00 %
objective value:                                    0
y                                                   2   (obj:0)
z                                     1.0000000000001   (obj:0)
alexhoen commented 1 year ago

First I would reply to the second point:

SimpleProbing treats z as constant resulting in the equation y = 2 - 2x and uses this equation to replace all appearances of y. This results in the constraint (2) to be -x >= -1 which is always true and since x has no further restrictions and no appearances in the objective, PaPILO chooses a "random" value for x. Afaik, this is valid and it also passes the SCIP feasibility check, so I don't see why x should be 1?

To 1: Epsilon is used to determine which values are zero. so if ub -lb <= epsilon then this means ub == lb and hence the variable should be fixed. Surely, this can result in some errors. You can always decrease epsilon to prevent this behavior or scale your matrix a priori to prevent this from happening.

liangbug commented 1 year ago

The true maximal object value of the simple model is 1, but gets 0. If you comment parameter setting code of the simple model, you will get $x=1$.

To 2.

Presolver simpleprobing applying 
row -10 col 1 val -2.0
row -1 col 0 val 2.0
tsx
modified rows: 0,1,
row -10 col 2 val -1.999511667349907e-13
row -1 col 0 val 1.0000000000001

According to the log, in ProblemUpdate.hpp ProblemUpdate::applyTransaction(), ColReduction::REPLACE, constraintMatrix.aggregate():

After that, in ProblemUpdate.hpp, ProblemUpdate::removeSingletonRow():

liangbug commented 1 year ago

To 1: I found the priority of Papilo presolvers in SCIP PresolMILP.cpp.

Fast    colsingleton
Fast    coefftightening
Fast    propagation
Medium  simpleprobing
Medium  parallelrows
Medium  stuffing
Medium  dualfix
Medium  fixcontinuous
Medium  simplifyineq
Medium  doubletoneq
Exhaustive  implint
Exhaustive  dualinfer
Exhaustive  probing
Exhaustive  domcol
Exhaustive  substitution

Because the priority of fixcontinuous is less than simpleprobing, the simpleprobing may get variables with 0 < ub-lb < epislon during several rounds of presolving. In my original problem, I get a variable with 0 < ub-lb < epislon in round 9 after probing, and then execute simpleprobing before fixcontinuous in round 11.

alexhoen commented 1 year ago

ah sorry, I am not that used to PySCIPOpt and I missed the objective function. Then there is definitely a problem there.

FYI: the priority of the presolvers can also be seen in the log where the results are displayed.

So, yes we should prioritize fixing continuous variables.

But, regardless, this should also return the correct results if disabled. I will take a look into it.

Thanks for your deep dive into PaPILO and SCIP. We really appreciate it.

alexhoen commented 1 year ago

fyi: a workaround would be to set the threads in PaPILO to more than 1 (presolving/milp/threads)

alexhoen commented 1 year ago

Could you please verify that with the newest master version of PaPILO SCIP returns the correct solution?

Thank you very much!!!

https://github.com/scipopt/papilo/commit/ebf0206ec42fa004f973ab9be5c17bb82d6de5ba

liangbug commented 1 year ago

I have verified that SCIP 8.0.3 with newest master version of PaPILO returns the correct solution. Even though the priority of simple probing is still higher than fixcontinus in SCIP presol_milp.cpp, but it doesn't matter. https://github.com/scipopt/scip/blob/5ecb558e560388c9786affaf5c9ca662c36cfd02/src/scip/presol_milp.cpp#L367 https://github.com/scipopt/scip/blob/5ecb558e560388c9786affaf5c9ca662c36cfd02/src/scip/presol_milp.cpp#L375

Thank you very much!!!

alexhoen commented 1 year ago

Perfect!! The bugfix will go live with the next bugfix release.

SCIP loads the presolving techniques and defines its own order (see https://www.scipopt.org/doc/html/presol__milp_8cpp_source.php line 346ff). Therefore, the order should stay and obviously stay the same. This was "only" the check if the presolving techniques simple-probing generates valid reductions. Regardless of the order of presolvers, this should be the case for all presolvers. We will change the order in SCIP probably with the next bugfix release

I will close this ticket now. If you have questions feel free to contact us!

Many thanks you too.