Conic duals: Primal infeasibilty certificate does not pass farkas test

conicdualfail.jl

using MAT, JLD
using Mosek
using MathProgBase

TOL = 1e-6
solver = MosekSolver()
case = load("conicduals.jld")
A = case["A"]
c = case["c"]
b = case["b"]
constr_cones = case["constr_cones"]
var_cones = case["var_cones"]

# constr_cone and the dual cone - K₁ & K₁*    
# 23 + 23 + 2 - :Zero cones -- Dual cones - :Free
# 37 + 37 - :NonNeg cones -- Dual cones - :NonNeg
# 37 :SOC (R×R^{10}) -- Dual cones - :SOC

# var_cones and the duals - K₂ & K₂*
# 48 :Free cones -- Dual cones - :Zero
# 38 :NonNeg cones -- Dual cones - :NonNeg

#= Primal and dual problems
(P) = Min c'x s.t.
b - Ax ∈ K₁
x ∈ K₂

(D) = Max - b'y s.t.
c + A'y ∈ K₂*
y ∈ K₁*
=#

#= Farkas certificate for primal infeasibility
A'y ∈ K₂*
y ∈ K₁*
- b'y > 0
=#

@assert length(c) == 48 + 38
@assert length(b) == 23 + 23 + 2 + 37 + 37 + 37*10

m = MathProgBase.model(solver)
MathProgBase.loadconicproblem!(m, c, A, b, constr_cones, var_cones)
MathProgBase.optimize!(m)
@show MathProgBase.status(m)
@show MathProgBase.getobjval(m)
y = MathProgBase.getconicdual(m)

# Farkas certificate check
@assert -dot(b,y) > 0
# Check if A'y ∈ K₂*
k2check = A'y
@assert all(k2check[1:48] .< TOL) && all(k2check[1:48] .> -TOL) 
@assert all(k2check[49:end] .>= -TOL)
# y ∈ K₁*
# first 48 elements are in the :Free cone - nothing to check
@assert all(y[49:48+37+37] .>= -TOL) # :NonNeg cone check
for i in 1:37
    index = 48 + 2*37 + (i-1)*(10)
    println(">> ||.|| = $(norm(y[index+2:index+10])), rhs = $(y[index+1])")
    @assert norm(y[index+2:index+10]) - y[index+1] <= TOL
end