No dual solution available for constraint in LP with presolve disabled.

[mjstahlberg]

I'm having issues with LP duals (with presolve disabled). Some constraints yield the correct duals, other constraints in the same problem and during the same optimization pass throw no dual solution available for constraint warnings. Setting the affected dual values to 0.0 is not correct; some of the affected constraints are supposed to have non-zero duals.

Unfortunately my example where some duals of tight constraints are returned while some are not is too involved for a bug report. However, here is an example with a single constraint where retrieving its dual fails even though the constraint is tight (and its dual value should be 1.0):

#!/usr/bin/python
import pyscipopt
m = pyscipopt.Model()
x = m.addVar("x")
c = m.addCons(x >= 1)
m.setObjective(x, "minimize")
m.setPresolve(pyscipopt.scip.PY_SCIP_PARAMSETTING.OFF)
m.optimize()
m.getDualsolLinear(c)

Output on v1.2.0:

Traceback (most recent call last):
 File "src/pyscipopt/scip.pyx", line 1276, in pyscipopt.scip.Model.getDualsolLinear
 File "src/pyscipopt/scip.pyx", line 1228, in pyscipopt.scip.Model.getTransformedCons
 File "src/pyscipopt/scip.pyx", line 283, in pyscipopt.scip.Constraint.create
Warning: cannot create Constraint with SCIP_CONS* == NULL

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
 File "./test_scip_duals.py", line 9, in <module>
 m.getDualsolLinear(c)
 File "src/pyscipopt/scip.pyx", line 1292, in pyscipopt.scip.Model.getDualsolLinear
Warning: no dual solution available for constraint c1

Output on current f5956e0:

Traceback (most recent call last):
 File "src/pyscipopt/scip.pyx", line 1577, in pyscipopt.scip.Model.getDualsolLinear
 File "src/pyscipopt/scip.pyx", line 1529, in pyscipopt.scip.Model.getTransformedCons
 File "src/pyscipopt/scip.pyx", line 405, in pyscipopt.scip.Constraint.create
Warning: cannot create Constraint with SCIP_CONS* == NULL

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
 File "./test_scip_duals.py", line 9, in <module>
 m.getDualsolLinear(c)
 File "src/pyscipopt/scip.pyx", line 1592, in pyscipopt.scip.Model.getDualsolLinear
Warning: no dual solution available for constraint c1

Python is 3.6.4, scipoptsuite is 5.0.1.

[mattmilten]

Could you please try disabling heuristics and propagation as well:

model.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
model.disablePropagation()

We had a similar problem where SCIP modified the problem internally and then the LP solver couldn't provide the dual solution anymore.

On your example model, it is already enough to disable propagation to get the correct result.

[mjstahlberg]

Hello @mattmilten, I can confirm that this fixes the problem with my small example. Unfortunately, just like with the fix that was proposed in #137, I am now getting -0.0 for the previously affected constraints in my more complex example. I'll try to get a plaintext version of it or create another one.

For disablePropagation, what would be the correct way to reset this to its default?

[mattmilten]

These are SCIP's default values for propagation:

self.setIntParam("propagating/maxroundsroot", 1000)
self.setIntParam("propagating/maxrounds", 100)

Can't you call writeProblem() before solving?

[mjstahlberg]

Here's my other example:

#!/usr/bin/python

import pyscipopt

M = pyscipopt.Model()

x = []
for i in range(5):
 x.append(M.addVar(name = "x[{}]".format(i), lb = -1e20, ub = 1e20))

M.setObjective( 1.0 * x[0] + 2.0 * x[1] + 3.0 * x[2] + 4.0 * x[3] + 5.0 * x[4], "maximize")

C = []
C.append(M.addCons(-1.0 * x[0] <= 1))
C.append(M.addCons( 1.0 * x[0] <= 1))
C.append(M.addCons(-2.0 * x[2] <= 1))
C.append(M.addCons( 2.0 * x[2] <= 1))
C.append(M.addCons(-2.0 * x[3] <= 1))
C.append(M.addCons( 2.0 * x[3] <= 1))
C.append(M.addCons(-1.0 * x[3] <= 1))
C.append(M.addCons( 1.0 * x[0] <= 1))
C.append(M.addCons( 1.0 * x[0] + 2.0 * x[2] <= 1))
C.append(M.addCons(-1.0 * x[0] - 2.0 * x[2] <= 1))
C.append(M.addCons( 1.0 * x[1] + 1.0 * x[2] + 1.0 * x[3] <= 1))
C.append(M.addCons(-1.0 * x[1] - 1.0 * x[2] - 1.0 * x[3] <= 1))
C.append(M.addCons( 1.0 * x[0] + 2.0 * x[1] <= 1))
C.append(M.addCons(-1.0 * x[0] - 2.0 * x[1] <= 1))
C.append(M.addCons( 1.0 * x[0] + 3.0 * x[2] + 1.0 * x[4] <= 1))
C.append(M.addCons(-1.0 * x[0] - 3.0 * x[2] - 1.0 * x[4] <= 1))

M.setPresolve(pyscipopt.SCIP_PARAMSETTING.OFF)
M.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
M.disablePropagation()
M.optimize()

print("Obj:", M.getObjVal())
for j in range(len(C)):
 print("C[{}]:".format(j), M.getDualsolLinear(C[j]))

The output is:

Obj: 14.0
C[0]: -0.0
C[1]: -0.0
C[2]: -0.0
C[3]: -0.0
C[4]: -0.0
C[5]: -0.0
C[6]: -0.0
C[7]: -0.0
C[8]: 0.0
C[9]: 5.0
C[10]: 0.0
C[11]: 0.0
C[12]: 1.0
C[13]: 0.0
C[14]: 5.0
C[15]: 0.0

Solving the same problem with Gurobi I get the same duals for C[8] to C[15] while C[3] = 1 and C[4] = 2 (all other duals are "positive" 0.0). Note that when you take adjacent pairs of constraints and sum up the absolute values of the difference of their duals you then get 14, the dual optimal value.

For completeness, here is the original test case written for PICOS. Both primal and dual tests pass for all of CPLEX, CVXOPT, GLPK, Gurobi, MOSEK and SMCP while for my work in progress reimplementation of SCIP only the primal test passes (the old SCIP implementation did not retrieve any duals):

 def setUp(self):
 # Data.
 c = picos.new_param("c", cvxopt.matrix([1,2,3,4,5]))
 A = picos.new_param("A", [
 cvxopt.matrix([1,0,0,0,0]).T,
 cvxopt.matrix([0,0,2,0,0]).T,
 cvxopt.matrix([0,0,0,2,0]).T,
 cvxopt.matrix([1,0,0,0,0]).T,
 cvxopt.matrix([1,0,2,0,0]).T,
 cvxopt.matrix([0,1,1,1,0]).T,
 cvxopt.matrix([1,2,0,0,0]).T,
 cvxopt.matrix([1,0,3,0,1]).T
 ])
 self.n = n = len(A)
 self.m = m = c.size[0]

 # Primal problem.
 self.P = P = picos.Problem()
 self.u = u = P.add_variable("u", m)
 P.set_objective("max", c|u)
 self.C = P.add_list_of_constraints(
 [abs(A[i]*u) < 1 for i in range(n)], "i", "[n]")

 # Dual problem.
 self.D = D = picos.Problem()
 self.z = z = [D.add_variable("z[{}]".format(i)) for i in range(n)]
 self.mu = mu = D.add_variable("mu", n)
 D.set_objective("min", 1|mu)
 D.add_list_of_constraints(
 [abs(z[i]) < mu[i] for i in range(n)], "i", "[n]")
 D.add_constraint(
 picos.sum([A[i].T * z[i] for i in range(n)], "i", "[n]") == c)

 self.skipUnsupported(P, D)

 def testPrimalSolution(self):
 self.primalSolve(self.P)
 self.expectObjective(self.P, 14.0)

 def testDualSolution(self):
 self.dualSolve(self.P)
 for i in range(self.n):
 self.z[i].set_value(self.C[i].dual[1] - self.C[i].dual[0])
 self.mu.set_value(cvxopt.matrix([abs(z_i.value) for z_i in self.z]))
 self.expectObjective(self.D, 14.0)

[mattmilten]

Sorry, I somehow lost track of this. Is this issue still relevant?

[mjstahlberg]

Tested again on current c317270. I get the same infeasible dual solution as before. Unfortunately this is still the only reproducible example that I have; if you are having trouble debugging it I could write a couple more basic LP test cases for PICOS, which I should do eventually anyway, and see if another one breaks.

[mattmilten]

I see, a basic LP test case would be much appreciated.

[mjstahlberg]

Sorry for the delay. I did find a smaller example where duals returned are incorrect:

#!/usr/bin/python

import pyscipopt

M = pyscipopt.Model()
x = M.addVar(name = "x", lb = -1e20, ub = 1e20)
M.setObjective(x, "maximize")
C1 = M.addCons(-1.0 * x <= 1)
C2 = M.addCons( 1.0 * x <= 1)

M.setPresolve(pyscipopt.SCIP_PARAMSETTING.OFF)
M.setHeuristics(pyscipopt.SCIP_PARAMSETTING.OFF)
M.disablePropagation()
M.optimize()

print()
print("Obj :", M.getObjVal())
print("x :", M.getVal(x))
print("Duals:", M.getDualsolLinear(C1), M.getDualsolLinear(C2))

The output is:

presolving:
presolving (0 rounds: 0 fast, 0 medium, 0 exhaustive):
 0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 0 cliques
presolved problem has 1 variables (0 bin, 0 int, 0 impl, 1 cont) and 2 constraints
 2 constraints of type <linear>
Presolving Time: 0.00

 time | node | left |LP iter|LP it/n| mem |mdpt |frac |vars |cons |cols |rows |cuts |confs|strbr| dualbound | primalbound | gap 
* 0.0s| 1 | 0 | 0 | - | 543k| 0 | - | 1 | 2 | 1 | 0 | 0 | 0 | 0 | 1.000000e+00 | 1.000000e+00 | 0.00%
 0.0s| 1 | 0 | 0 | - | 543k| 0 | - | 1 | 2 | 1 | 0 | 0 | 0 | 0 | 1.000000e+00 | 1.000000e+00 | 0.00%

SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes : 1
Primal Bound : +1.00000000000000e+00 (1 solutions)
Dual Bound : +1.00000000000000e+00
Gap : 0.00 %

Obj : 1.0
x : 1.0
Duals: 1.0 1.0

SCIP claims that x = 1 is an optimum primal solution, which is clearly correct. However, complementary slackness does not hold as the first constraint -1.0 * x <= 1 is not tight but the dual returned for it is nonzero.

[mattmilten]

Thank you for providing such a small example for reproducing the error!

[fserra]

oh so the problem in our code was that we always assumed that in the constraints representing bounds the coefficient was one?