Branch-and-bound node incorrectly declared infeasible for bc1_presolved

[akazachk]

When running instance bc1_presolved, using option -d 64, I eventually get an error

    ## Solving for parent node 1/32. ##
    Clp0006I 0  Obj 1.7567017 Primal inf 14880.59 (44) Dual inf 383.62191 (184)
    Clp0003I Stopped - objective value 1.488059e+14
    Clp0003I Stopped - objective value 3.1825633e+14
    *** ERROR (version #eb19147): src/branch/VPCEventHandler.cpp:1281: Solver not proven optimal for node 90.

I try to replicate the partial tree generated by the vpc code using a smaller (easier to debug) code: bb_example.
Since the coin_example repository is currently private, I am also attaching the bb_example.zip code here.

The bb_example code proceeds identically to the vpc code in generating the partial branch-and-bound tree, up to node 30.
At that point, the underlying Clp calls diverge.

As shown in the screen shot, the issue is that on the left (the bb_example code), node with id 31 is declared optimal, while on the right (the vpc code), the same node is declared primal infeasible.

This node corresponds to setting the following bounds: (947,u,0), (963,u,0), (964,u,0), (965,l,1), (977,u,0).

As shown in the screenshot, Gurobi finds this node feasible, with an objective close to the one that Clp finds in the bb_example code.

---

Compilation of COIN code to reproduce:

1. get coinbrew script
2. install with

  chmod u+x coinbrew
  ./coinbrew fetch Cbc@master
  ./coinbrew build Cbc -b buildg -p buildg --no-prompt --enable-debug ADD_CXXFLAGS="-DSAVE_NODE_INFO" --with-cplex=false

[akazachk]

Adding a little bit more detail, some of the above is a printing artifact. It seems like Clp actually proceeds the same for much of the solution of the node with id 31 (in both the bb_example and vpc codes):

There is divergence at iteration 752 (possibly earlier, as I have not checked in more detail yet).

---

To generate more details for Clp printing within Cbc, you can set:

p ((OsiClpSolverInterface*)cbc_model->solver())->getModelPtr()->handler_->logLevel_=3
p ((OsiClpSolverInterface*)cbc_model->solver())->getModelPtr()->minIntervalProgressUpdate_=-1

[akazachk]

In the beginning, we are using ClpSimplexDual, called repeatedly from ClpSimplex.cpp:5503, and at some point we switch to ClpSimplexPrimal. I think this happens at line ClpSimplex.cpp:5623 when problemStatus_ = 10, at numberIterations_ = 304 (maybe earlier?). Eventually, the two codes diverge.

The first time I detect a difference is when numberIterations_ is 460. At line ClpSimplexPrimal.cpp:999, before gutsOfSolution is called, both codes have sumDualInfeasibilities_ = 2.6513519195572522e+17. However, in this iteration, the solution_ vectors start to diverge subtly (but why?). We get that sumDualInfeasibilities_ = 54960487247979456 in the vpc code, and it equals 54960487247938216 in the bb_example code.

Eventually, at this same line, we get to a substantially different solution in the two codes, if we set a breakpoint for when numberIterations_ is at least 752. Continuing checking this breakpoint, somewhere in an iteration past 1450, sumDualInfeasibilties_ goes to 0 in the bb_example code, but it might not in the vpc code. Also solution_[1] grows very large in the vpc code (~10^12), while it remains a reasonable value (in the 100s) for the bb_example code, growing smaller until it eventually is 0, around numberIterations_ >= 1680, and afterwards an optimal objective is returned for the bb_example_code, while the vpc one continues to struggle for another couple hundred iterations.

Why does this happen?

---

Debugging steps (in gdb) for the current versions of the two codes

vpc:

br VPCEventHandler.cpp:455 if numLeafNodes_ >= 22
run -f bc1_presolved.mps.gz -d 64
  [wait for breakpoint]
br ClpSimplexPrimal.cpp:999 if numberIterations >= 460
  [wait for breakpoint]
  999         numberThrownOut = gutsOfSolution(NULL, NULL, sayValuesPass);
p sumDualInfeasibilities_
  54960487247979456
p solution_[1]
  -6947.9237437694092
cond 1 numberIterations_ >= 752
  [wait for breakpoint]
 p nonLinearCost_->feasibleReportCost()
  73.285936541638094
next
  1000        double sumInfeasibility = nonLinearCost_->sumInfeasibilities();
call (int) printf("obj: %.20g\nsumDualInfeasibilities_ = %.20g\nsolution_[1] = %f\n", nonLinearCost_->feasibleReportCost(), sumDualInfeasibilities_, solution_[1])
  obj: 76.41460699991704075273
  sumDualInfeasibilities_ = 7785816558235727
  solution_[1] = -8215.733952

bb_example:

br bb_example.cpp:79 if numLeafNodes_ >= 22
run bc1_presolved.mps.gz
  [wait for breakpoint]
br ClpSimplexPrimal.cpp:999 if numberIterations >= 460
  [wait for breakpoint]
  999         numberThrownOut = gutsOfSolution(NULL, NULL, sayValuesPass);
p sumDualInfeasibilities_
  54960487247938216
p solution_[1]
  -6947.923743769571
cond 1 numberIterations_ >= 752
  [wait for breakpoint]
 p nonLinearCost_->feasibleReportCost()
  73.285936541638094
next
  1000        double sumInfeasibility = nonLinearCost_->sumInfeasibilities();
call (int) printf("obj: %.20g\nsumDualInfeasibilities_ = %.20g\nsolution_[1] = %f\n", nonLinearCost_->feasibleReportCost(), sumDualInfeasibilities_, solution_[1])
  obj: 53.51715626364362776712
  sumDualInfeasibilities_ = 319579030469625024
  solution_[1] = -6398.854461