Detect primal/dual infeasibility

src/batch/madipm/kernels.jl

@@ -14,6 +14,84 @@ function dual_objective!(dual_obj, y_vals, rhs_vals, zl_r, xl_r, zu_r, xu_r,
     return dual_obj
 end
 
+function update_primal_infeasibility_certificate!(batch_solver::AbstractBatchMPCSolver{T}) where T
+    ws = batch_solver.workspace
+    y = MadNLP.full(batch_solver.y)
+    rhs = MadNLP.full(batch_solver.rhs)
+    jacl = MadNLP.full(batch_solver.jacl)
+    zl = MadNLP.full(batch_solver.zl)
+    zu = MadNLP.full(batch_solver.zu)
+    nlb, nub = batch_solver.d.nlb, batch_solver.d.nub
+
+    cert_norm = ws.sum_lb
+    batch_mapreduce!(abs, max, zero(T), cert_norm, y)
+    if nlb > 0
+        batch_mapreduce!(abs, max, zero(T), ws.sum_ub, lower(batch_solver.zl))
+        @. cert_norm = max(cert_norm, ws.sum_ub)
+    end
+    if nub > 0
+        batch_mapreduce!(abs, max, zero(T), ws.sum_ub, upper(batch_solver.zu))
+        @. cert_norm = max(cert_norm, ws.sum_ub)
+    end
+    @. cert_norm = max(one(T), cert_norm)
+
+    batch_mapreduce!((jl, zl_i, zu_i) -> abs(jl - zl_i + zu_i), max, zero(T),
+                     ws.primal_cert_res, jacl, zl, zu)
+    @. ws.primal_cert_res /= cert_norm
+
+    data_norm = ws.sum_ub
+    batch_mapreduce!(abs, max, zero(T), data_norm, rhs)
+    if nlb > 0
+        batch_mapreduce!(abs, max, zero(T), ws.mu_curr, lower(batch_solver.xl))
+        @. data_norm = max(data_norm, ws.mu_curr)
+    end
+    if nub > 0
+        batch_mapreduce!(abs, max, zero(T), ws.mu_curr, upper(batch_solver.xu))
+        @. data_norm = max(data_norm, ws.mu_curr)
+    end
+    @. data_norm = max(one(T), data_norm)
+    @. ws.primal_cert_margin = ws.dual_obj / (cert_norm * data_norm)
+    return
+end
+
+function update_dual_infeasibility_certificate!(batch_solver::AbstractBatchMPCSolver{T}) where T
+    ws = batch_solver.workspace
+    if !batch_solver.bcb.nlp.meta.islp
+        fill!(ws.dual_cert_res, typemax(T))
+        fill!(ws.dual_cert_bound, typemax(T))
+        fill!(ws.dual_cert_margin, -typemax(T))
+        return
+    end
+
+    x = MadNLP.primal(batch_solver.x)
+    ray = MadNLP.primal(batch_solver._w1)
+    batch_mapreduce!(abs, max, zero(T), ws.sum_lb, x)
+    @. ws.sum_lb = max(one(T), ws.sum_lb)
+    @. ray = x / ws.sum_lb
+
+    batch_mapreduce!((c_i, rhs_i) -> abs(c_i + rhs_i), max, zero(T),
+                     ws.dual_cert_res, MadNLP.full(batch_solver.c), MadNLP.full(batch_solver.rhs))
+    @. ws.dual_cert_res /= ws.sum_lb
+
+    nlb, nub = batch_solver.d.nlb, batch_solver.d.nub
+    fill!(ws.dual_cert_bound, zero(T))
+    if nlb > 0
+        batch_mapreduce!(d -> max(zero(d), -d), max, zero(T),
+                         ws.dual_cert_bound, xp_lr(batch_solver._w1))
+    end
+    if nub > 0
+        batch_mapreduce!(d -> max(zero(d), d), max, zero(T),
+                         ws.sum_ub, xp_ur(batch_solver._w1))
+        @. ws.dual_cert_bound = max(ws.dual_cert_bound, ws.sum_ub)
+    end
+
+    batch_mapreduce!(*, +, zero(T), ws.dual_cert_margin, MadNLP.primal(batch_solver.f), ray)
+    batch_mapreduce!(abs, max, zero(T), ws.sum_ub, MadNLP.primal(batch_solver.f))
+    @. ws.sum_ub = max(one(T), ws.sum_ub)
+    @. ws.dual_cert_margin = -ws.dual_cert_margin / ws.sum_ub
+    return
+end
+
 function set_initial_primal_rhs!(solver::AbstractBatchMPCSolver)
     p = solver.p
     fill!(MadNLP.full(p), 0.0)

src/batch/madipm/solver.jl

@@ -245,6 +245,11 @@ function initialize_solver_state!(batch_solver::AbstractBatchMPCSolver{T}) where
     fill!(ws.inf_du, zero(T))
     fill!(ws.inf_compl, zero(T))
     fill!(ws.dual_obj, zero(T))
+    fill!(ws.primal_cert_res, typemax(T))
+    fill!(ws.primal_cert_margin, -typemax(T))
+    fill!(ws.dual_cert_res, typemax(T))
+    fill!(ws.dual_cert_bound, typemax(T))
+    fill!(ws.dual_cert_margin, -typemax(T))
     fill!(ws.alpha_p, zero(T))
     fill!(ws.alpha_d, zero(T))
     t_now = time()
@@ -267,20 +272,36 @@ function compute_term_gpu!(ws::UniformBatchWorkspace{T}, opt) where T
     Int_INFEASIBLE = Int(MadNLP.INFEASIBLE_PROBLEM_DETECTED)
     Int_DIVERGING = Int(MadNLP.DIVERGING_ITERATES)
     Int_REGULAR = Int(MadNLP.REGULAR)
+    cert_enabled = opt.certificate_termination
+    pcert_tol = T(opt.primal_infeasibility_cert_tol)
+    dcert_tol = T(opt.dual_infeasibility_cert_tol)
     @. ws._term_gpu = ifelse(
         ws._ls_error > zero(Int32),
         Int_ERROR,
         ifelse(
             max(ws.inf_pr, ws.inf_du, ws.inf_compl) <= tol,
             Int_SOLVED,
             ifelse(
-                (ws.inf_compl > div_tol * ws.best_complementarity) &
-                (ws.dual_obj > max(ds * abs(ws.obj_val), one(T))),
+                cert_enabled &
+                (ws.primal_cert_res <= pcert_tol) &
+                (ws.primal_cert_margin > pcert_tol),
                 Int_INFEASIBLE,
                 ifelse(
-                    ws.obj_val < -(div_tol * max(ds * abs(ws.dual_obj), one(T))),
+                    cert_enabled &
+                    (ws.dual_cert_res <= dcert_tol) &
+                    (ws.dual_cert_bound <= dcert_tol) &
+                    (ws.dual_cert_margin > dcert_tol),
                     Int_DIVERGING,
-                    Int_REGULAR,
+                    ifelse(
+                        (ws.inf_compl > div_tol * ws.best_complementarity) &
+                        (ws.dual_obj > max(ds * abs(ws.obj_val), one(T))),
+                        Int_INFEASIBLE,
+                        ifelse(
+                            ws.obj_val < -(div_tol * max(ds * abs(ws.dual_obj), one(T))),
+                            Int_DIVERGING,
+                            Int_REGULAR,
+                        ),
+                    ),
                 ),
             ),
         ),
@@ -313,6 +334,8 @@ function update_termination_criteria!(batch_solver::AbstractBatchMPCSolver{T}) w
         lower(zl), lower(xl), upper(zu), upper(xu),
         ws.sum_lb, ws.sum_ub, nlb, nub)
 
+    update_primal_infeasibility_certificate!(batch_solver)
+    update_dual_infeasibility_certificate!(batch_solver)
     compute_term_gpu!(ws, opt)
     return
 end

src/batch/structure.jl

@@ -21,6 +21,11 @@ struct UniformBatchWorkspace{T, VT<:AbstractVector{T}, MT<:AbstractMatrix{T}, MI
     inf_compl::MT
     best_complementarity::MT
     dual_obj::MT
+    primal_cert_res::MT
+    primal_cert_margin::MT
+    dual_cert_res::MT
+    dual_cert_bound::MT
+    dual_cert_margin::MT
     status::Vector{MadNLP.Status}
 
     _term_gpu::MI64
@@ -61,6 +66,11 @@ function UniformBatchWorkspace(::Type{MT}, ::Type{VT}, n::Int, m::Int, nlb::Int,
         MT(undef, 1, batch_size),  # inf_compl
         MT(undef, 1, batch_size),  # best_complementarity
         MT(undef, 1, batch_size),  # dual_obj
+        MT(undef, 1, batch_size),  # primal_cert_res
+        MT(undef, 1, batch_size),  # primal_cert_margin
+        MT(undef, 1, batch_size),  # dual_cert_res
+        MT(undef, 1, batch_size),  # dual_cert_bound
+        MT(undef, 1, batch_size),  # dual_cert_margin
         fill(MadNLP.INITIAL, batch_size),  # status
         similar(_proto, Int64, 1, batch_size),    # _term_gpu
         zeros(Int64, 1, batch_size),  # _term_cpu

src/kernels.jl

@@ -437,6 +437,58 @@ function dual_objective(solver::MPCSolver)
     return dobj
 end
 
+function has_primal_infeasibility_certificate(solver::MPCSolver{T}) where {T}
+    solver.opt.certificate_termination || return false
+    cert_norm = max(
+        one(T),
+        norm(solver.y, Inf),
+        norm(solver.zl_r, Inf),
+        norm(solver.zu_r, Inf),
+    )
+    data_norm = max(
+        one(T),
+        norm(solver.rhs, Inf),
+        norm(solver.xl_r, Inf),
+        norm(solver.xu_r, Inf),
+    )
+    residual = MadNLP.primal(solver._w2)
+    copyto!(residual, solver.jacl)
+    residual .-= MadNLP.full(solver.zl)
+    residual .+= MadNLP.full(solver.zu)
+    cert_res = norm(residual, Inf) / cert_norm
+    cert_margin = dual_objective(solver) / (cert_norm * data_norm)
+    tol = solver.opt.primal_infeasibility_cert_tol
+    return isfinite(cert_res) && isfinite(cert_margin) &&
+           cert_res <= tol && cert_margin > tol
+end
+
+function has_dual_infeasibility_certificate(solver::MPCSolver{T}) where {T}
+    solver.opt.certificate_termination || return false
+    solver.class isa LinearProgram || return false
+    # For QPs, we also need to check the curvative, but this is not implemented.
+    # So we only detect dual infeasibility for LPs. For QPs, we rely on the
+    # diverging iterates check.
+
+    x = MadNLP.primal(solver.x)
+    ray_norm = max(one(T), norm(x, Inf))
+
+    lhs = MadNLP.dual(solver._w2)
+    copyto!(lhs, solver.c)
+    lhs .+= solver.rhs
+    cert_res = norm(lhs, Inf) / ray_norm
+
+    lower_violation = isempty(solver.x_lr) ? zero(T) : max(zero(T), -minimum(solver.x_lr) / ray_norm)
+    upper_violation = isempty(solver.x_ur) ? zero(T) : max(zero(T), maximum(solver.x_ur) / ray_norm)
+    bound_violation = max(lower_violation, upper_violation)
+
+    obj_ray = dot(MadNLP.primal(solver.f), x) / ray_norm
+    obj_norm = max(one(T), norm(MadNLP.primal(solver.f), Inf))
+    cert_margin = -obj_ray / obj_norm
+    tol = solver.opt.dual_infeasibility_cert_tol
+    return isfinite(cert_res) && isfinite(bound_violation) && isfinite(cert_margin) &&
+           cert_res <= tol && bound_violation <= tol && cert_margin > tol
+end
+
 function get_optimality_gap(solver::MPCSolver)
     return MadNLP.get_inf_compl(
         solver.x_lr,

src/solver.jl

@@ -206,6 +206,10 @@ function update_termination_criteria!(solver::MadNLP.AbstractMadNLPSolver)
 
     if max(solver.inf_pr, solver.inf_du, solver.inf_compl) <= solver.opt.tol
         solver.status = MadNLP.SOLVE_SUCCEEDED
+    elseif has_primal_infeasibility_certificate(solver)
+        solver.status = MadNLP.INFEASIBLE_PROBLEM_DETECTED
+    elseif has_dual_infeasibility_certificate(solver)
+        solver.status = MadNLP.DIVERGING_ITERATES  # TODO: MadNLP.UNBOUNDED_PROBLEM_DETECTED?
     elseif ((solver.inf_compl > solver.opt.divergence_tol * solver.best_complementarity) &&
             (dobj > max(solver.opt.divergence_scale * abs(solver.obj_val), 1.0)))
         solver.status = MadNLP.INFEASIBLE_PROBLEM_DETECTED

src/structure.jl

@@ -76,6 +76,24 @@ mutable struct MPCSolver{
     status::MadNLP.Status
 end
 
+for (k, attribute) in enumerate(fieldnames(MPCSolver))
+    fname = "get_$(attribute)"
+    sf = Symbol(fname)
+    if isdefined(MadNLP, sf)
+        @eval begin
+            @inline function MadNLP.$(sf)(solver::MPCSolver)
+                return getfield(solver, $k)
+            end
+        end
+    else
+        @eval begin
+            @inline function $(sf)(solver::MPCSolver)
+                return getfield(solver, $k)
+            end
+        end
+    end
+end
+
 function MPCSolver(nlp::NLPModels.AbstractNLPModel{T,VT}; kwargs...) where {T, VT}
     options = load_options(nlp; kwargs...)
 
@@ -185,4 +203,3 @@ function MadNLP.print_iter(solver::MPCSolver; options...)
         solver.alpha_d,solver.alpha_p))
     return
 end
-

src/utils.jl

@@ -80,6 +80,9 @@ end
     max_wall_time::Float64 = 1e6
     divergence_tol::Float64 = 1e4
     divergence_scale::Float64 = 10.0
+    certificate_termination::Bool = true
+    primal_infeasibility_cert_tol::Float64 = 1e-7
+    dual_infeasibility_cert_tol::Float64 = 1e-7
     kappa_d::Float64 = 1e-5
     fixed_variable_treatment::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? MadNLP.RelaxBound : MadNLP.MakeParameter
     equality_treatment::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? MadNLP.RelaxEquality : MadNLP.EnforceEquality

test/batch/solver.jl

@@ -331,6 +331,11 @@ end
         ws.best_complementarity .= best
         ws.obj_val .= obj
         ws.dual_obj .= dobj
+        ws.primal_cert_res .= 1.0
+        ws.primal_cert_margin .= -1.0
+        ws.dual_cert_res .= 1.0
+        ws.dual_cert_bound .= 1.0
+        ws.dual_cert_margin .= -1.0
         fill!(ws._ls_error, ls)
     end
     _s(j) = MadNLP.Status(ws._term_gpu[1,j])
@@ -347,6 +352,19 @@ end
     MadIPM.compute_term_gpu!(ws, bat.opt)
     @test _s(1) == MadNLP.INFEASIBLE_PROBLEM_DETECTED
 
+    _set!()
+    ws.primal_cert_res .= bat.opt.primal_infeasibility_cert_tol / 10
+    ws.primal_cert_margin .= bat.opt.primal_infeasibility_cert_tol * 10
+    MadIPM.compute_term_gpu!(ws, bat.opt)
+    @test _s(1) == MadNLP.INFEASIBLE_PROBLEM_DETECTED
+
+    _set!()
+    ws.dual_cert_res .= bat.opt.dual_infeasibility_cert_tol / 10
+    ws.dual_cert_bound .= bat.opt.dual_infeasibility_cert_tol / 10
+    ws.dual_cert_margin .= bat.opt.dual_infeasibility_cert_tol * 10
+    MadIPM.compute_term_gpu!(ws, bat.opt)
+    @test _s(1) == MadNLP.DIVERGING_ITERATES
+
     _set!(dobj=1.0, obj=-(dt*ds*2))
     MadIPM.compute_term_gpu!(ws, bat.opt)
     @test _s(1) == MadNLP.DIVERGING_ITERATES

test/runtests.jl

@@ -156,6 +156,50 @@ end
         @test flag
     end
 
+    @testset "Certificate termination" begin
+        infeas_lp = QuadraticModel(
+            [0.0],
+            Int[],
+            Int[],
+            Float64[];
+            Arows = [1],
+            Acols = [1],
+            Avals = [1.0],
+            lcon = [-1.0],
+            ucon = [-1.0],
+            lvar = [0.0],
+            uvar = [Inf],
+            x0 = [0.0],
+        )
+        infeas_solver = MadIPM.MPCSolver(infeas_lp; print_level=MadNLP.ERROR, scaling=false)
+        @test infeas_solver.class isa MadIPM.LinearProgram
+        MadIPM.initialize!(infeas_solver)
+        infeas_solver.y .= 1.0
+        infeas_solver.zl_r .= 1.0
+        MadNLP.jtprod!(infeas_solver.jacl, infeas_solver.kkt, infeas_solver.y)
+        @test MadIPM.has_primal_infeasibility_certificate(infeas_solver)
+        MadIPM.update_termination_criteria!(infeas_solver)
+        @test infeas_solver.status == MadNLP.INFEASIBLE_PROBLEM_DETECTED
+
+        unbounded_lp = QuadraticModel(
+            [-1.0],
+            Int[],
+            Int[],
+            Float64[];
+            lvar = [0.0],
+            uvar = [Inf],
+            x0 = [1.0],
+        )
+        unbounded_solver = MadIPM.MPCSolver(unbounded_lp; print_level=MadNLP.ERROR, scaling=false)
+        @test unbounded_solver.class isa MadIPM.LinearProgram
+        MadIPM.initialize!(unbounded_solver)
+        MadNLP.primal(unbounded_solver.x) .= 1.0
+        MadIPM.evaluate_model!(unbounded_solver)
+        @test MadIPM.has_dual_infeasibility_certificate(unbounded_solver)
+        MadIPM.update_termination_criteria!(unbounded_solver)
+        @test unbounded_solver.status == MadNLP.DIVERGING_ITERATES
+    end
+
     @testset "Standard formulation" begin
         new_qp = MadIPM.standard_form_qp(qp)
         solver = MadIPM.MPCSolver(new_qp; print_level=MadNLP.ERROR)