From 9dfb91e5c7321dd14f445fd4e93cc8b8e23ca903 Mon Sep 17 00:00:00 2001 From: "Klamkin, Michael" Date: Wed, 20 May 2026 13:40:55 -0400 Subject: [PATCH 1/2] generate getters for madnlp compat --- src/structure.jl | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/structure.jl b/src/structure.jl index ed492057..bdbf6088 100644 --- a/src/structure.jl +++ b/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 - From a3bd24c6b2e3dd94990504fdd7fe36d86ea82974 Mon Sep 17 00:00:00 2001 From: "Klamkin, Michael" Date: Wed, 20 May 2026 14:03:22 -0400 Subject: [PATCH 2/2] add certificates --- src/batch/madipm/kernels.jl | 78 +++++++++++++++++++++++++++++++++++++ src/batch/madipm/solver.jl | 31 +++++++++++++-- src/batch/structure.jl | 10 +++++ src/kernels.jl | 52 +++++++++++++++++++++++++ src/solver.jl | 4 ++ src/utils.jl | 3 ++ test/batch/solver.jl | 18 +++++++++ test/runtests.jl | 44 +++++++++++++++++++++ 8 files changed, 236 insertions(+), 4 deletions(-) diff --git a/src/batch/madipm/kernels.jl b/src/batch/madipm/kernels.jl index 185163a9..55e96157 100644 --- a/src/batch/madipm/kernels.jl +++ b/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) diff --git a/src/batch/madipm/solver.jl b/src/batch/madipm/solver.jl index 840c047d..f8fa958d 100644 --- a/src/batch/madipm/solver.jl +++ b/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,6 +272,9 @@ 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, @@ -274,13 +282,26 @@ function compute_term_gpu!(ws::UniformBatchWorkspace{T}, opt) where T 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 diff --git a/src/batch/structure.jl b/src/batch/structure.jl index 44db0b73..bbeddaa6 100644 --- a/src/batch/structure.jl +++ b/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 diff --git a/src/kernels.jl b/src/kernels.jl index 88cd5687..07465685 100644 --- a/src/kernels.jl +++ b/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, diff --git a/src/solver.jl b/src/solver.jl index 8d8d32cc..4b53d1eb 100644 --- a/src/solver.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index d8773830..0c0de3e2 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/batch/solver.jl b/test/batch/solver.jl index 3f49f09e..a5449562 100644 --- a/test/batch/solver.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 3213ddb6..2fa23635 100644 --- a/test/runtests.jl +++ b/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)