Interpreting infeasibility/unboundedness ray in CLP

[febattista]

Hi, I'm developing functionalities in another COIN-OR project, namely the SYMPHONY MILP branch-and-cut solver, which is using CLP (via OSI) for solving LP relaxations.

I need to collect a ray proving primal infeasibility/dual unboundedness, whenever a LP is primal infeasible. As of now, I achieve that by calling the method OsiClpSolverInterface.getDualRays(). The issue is that, apparently, this function returns a non-null vector<double> only when the method OsiClpSolverInterface.initialSolve() is invoked for solving the LP, while for subsequent calls of OsiClpSolverInterface.resolve() getDualRays() will return an empty vector even if isProvenDualInfeasible() == true. It is important in my context that a ray is always returned for infeasible primal LPs. As of now, when this happens I'm resolving from scratch an additional LP by calling initialSolve(), but of course I'd like not to undertake this additional computational burden. Would it be possible to achieve that?

Given an LP in general form (i.e. with any combination of <=, >=, = constraints), defined by matrix $A$, right-hand side $b$ and objective func $c$, with general bounds $l \leq x \leq u$ given as input to CLP and suppose this LP is primal infeasible. Given the ray $\lambda$ returned by CLP, I need to verify the Farkas proof apply for a different set of $b^\prime, l^\prime, u^\prime$. My question is the following: Does the ray $\lambda$ returned by CLP correspond to a ray for the LP in general form as given in input? Or does it correspond to a ray to the LP CLP is internally creating (e.g. after preprocessing/scaling...)? Clearly, I would need the ray to correspond to the LP as given in input.

P.S. these are the option being set:

• setupForRepeatedUse()
• setHintParam(OsiDoScale,false,OsiHintDo)
• setHintParam(OsiDoPresolveInInitial, false, OsiHintDo)
• setHintParam(OsiDoPresolveInResolve, false, OsiHintDo) (only when resolve() is called)
• setCleanupScaling(1)
• int sp = si->specialOptions(); if ((sp & 2) != 0) sp ^= 2; si->setSpecialOptions(sp); to set special options.

Thanks for you time in advance, Federico

[jjhforrest]

The resolve issue is easy to fix. For some reason OsiClp thinks its default is in branch and bound. Resolve takes a simpler route than initialSolve which does create the ray. Standard dual thinks it is in branch and bound and so does not. However there is an option to say create a ray even in branch and bound. To activate this add

int options = si.getModelPtr()->specialOptions(); si.getModelPtr()->setSpecialOptions(options|32);

You may wish to do si.getModelPtr()->setSpecialOptions(options); later to restore default.

For the scaling and presolve issues, I may need examples where it is wrong.

[febattista]

Thank for your answer John.

Regarding the ray, here's the issue I'm encountering.

At some point during the B&B in SYMPHONY, CLP solves the following infeasible LP:

\Problem name: ClpDefaultName

Minimize obj: -895 x0 -629 x1 -105 x2 -976 x3 -672 x4 -508 x5 -16 x6 -216 x7 -235 x8 -192 x9 Subject To cons0: -431 x0 -548 x1 -16 x2 -39 x3 -25 x4 -922 x5 -521 x6 -648 x7 -92 x8 -522 x9 <= -2747 cons1: -541 x0 -131 x1 -57 x2 -964 x3 -385 x4 -866 x5 -379 x6 -662 x7 -508 x8 -562 x9 <= -2740.50000 cons2: 665 x0 + 402 x1 + 191 x2 + 769 x3 + 378 x4 + 862 x5 + 277 x6 + 394 x7 + 77 x8 + 512 x9 + x10 = 2263 Bounds 0 <= x0 <= 0 0 <= x1 <= 1 0 <= x2 <= 0 0 <= x3 <= 0 0 <= x4 <= 0 1 <= x5 <= 1 0 <= x6 <= 0 0 <= x7 <= 1 0 <= x8 <= 1 0 <= x9 <= 1 End

By calling si->getDualRays(1, false);, CLP returns $\lambda = (0, -0.03957..., 0)$. Since I want $\lambda$ to be normalized, I divide it by its norm. Hence $\lambda = (0, -1, 0)$. Then I compute myself $\lambda A$ by

• const CoinPackedMatrix *A = si->getMatrixByCol();
• A->transposeTimes(ray, rayA);

Here, $\lambda A = (541, 131, 57, 964, 385, 866, 379, 662, 508, 562, 0)$.

Then, at some point in my code I want to check if the Farkas proof applies to the following LP (spoiler: it is feasible):

\Problem name: ClpDefaultName

Minimize obj: -895 x0 -629 x1 -105 x2 -976 x3 -672 x4 -508 x5 -16 x6 -216 x7 -235 x8 -192 x9 Subject To cons0: -431 x0 -548 x1 -16 x2 -39 x3 -25 x4 -922 x5 -521 x6 -648 x7 -92 x8 -522 x9 <= -2060.25000 cons1: -541 x0 -131 x1 -57 x2 -964 x3 -385 x4 -866 x5 -379 x6 -662 x7 -508 x8 -562 x9 <= -2349 cons2: 665 x0 + 402 x1 + 191 x2 + 769 x3 + 378 x4 + 862 x5 + 277 x6 + 394 x7 + 77 x8 + 512 x9 + x10 = 2263 Bounds 1 <= x0 <= 1 0 <= x1 <= 0 0 <= x2 <= 1 0 <= x3 <= 0 1 <= x4 <= 1 0 <= x5 <= 0 1 <= x6 <= 1 1 <= x7 <= 1 0 <= x8 <= 0 1 <= x9 <= 1 End

The Farkas proof I'm testing goes as $\lambda A z - \lambda b > 0$, where $z_j = l_j$, if $\lambda A_j \geq 0$ and $z_j = u_j$, if $\lambda A_j < 0$, for $j=0,...,n-1$, yield the minimum violation of the Farkas constraint. Here since $\lambda A \geq 0$, we have $z = l$.

In this case, we have: $\lambda A z = 541 + 385 + 379 + 662 + 562 = 2529$ and $\lambda b = 2349$, which yields $\lambda A z - \lambda b = 2529 - 2349 = 180 > 0$. So, according to this Farkas proof, this LP should be infeasible, but solving it with CLP one can see that this is not the case.

[febattista]

Of course, I assume the ray from CLP is correct and I'm more likely to think that my way of testing Farkas proof is wrong (probably just a matter of signs). This is why I asked to which form of LP the ray from CLP corresponds to or what form of Farkas proof I should apply. That proof test came from GUROBI's reference manual https://www.gurobi.com/documentation/current/refman/farkasproof.html which seems to be exactly what I'm looking for this LP general form.

The proof works correctly for most cases in my code:

1. I compute the Farkas proof as stated before (say farkasProof = true | false),
2. I solve the corresponding LP in CLP from scratch (of course I don't want to do it when I figure out how to make this work),
3. if farkasProof == true, I assert si->isProvenPrimalInfeasible() In some cases the assertion fails, as in the example above.