From 6e68858d513a1171e92c7e4f16f614a36d5183d9 Mon Sep 17 00:00:00 2001 From: GhislainDv Date: Tue, 25 Aug 2026 16:19:24 +0200 Subject: [PATCH 1/2] first draft for active set detection --- Project.toml | 2 + src/MadCheck.jl | 1 + src/active_set.jl | 693 +++++++++++++++++++++++++++++++++++++++++++++- test/Project.toml | 1 + test/runtests.jl | 24 ++ 5 files changed, 720 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7742fb0..23710d3 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,14 @@ version = "0.1.0" authors = ["fpacaud "] [deps] +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] +JuMP = "1.31.2" LinearAlgebra = "1.12.0" MadNLP = "0.10.1" NLPModels = "0.21.12" diff --git a/src/MadCheck.jl b/src/MadCheck.jl index 065fa2d..ebe3c7e 100644 --- a/src/MadCheck.jl +++ b/src/MadCheck.jl @@ -3,6 +3,7 @@ module MadCheck using LinearAlgebra using SparseArrays import NLPModels +using JuMP include("solution.jl") include("constraint_qualification.jl") diff --git a/src/active_set.jl b/src/active_set.jl index 2c37f56..9f477e6 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -1,4 +1,695 @@ +""" + AbstractActiveSetMethod -# TODO: return current estimation of the active set based on Oberlin & Wright 2005 +Base type for the different active set methods implemented +""" +abstract type AbstractActiveSetMethod end +""" + ActiveMethodSimple <: AbstractActiveSetMethod + +composite type for the naive method consisting of runing the test c(x) = 0, +contains one parameter + `tol`: tolerance of the method +""" + +struct ActiveMethodSimple <: AbstractActiveSetMethod + tol::Float64 +end + +""" + ActiveMethodSimple(; kwargs...) + +Creates an ActiveMethodeSimple with `tol` = 1e-8. + +""" + +ActiveMethodSimple(;tol = 1e-8) = ActiveMethodSimple(tol) + +""" + ActiveMethodLP_P <: AbstractActiveSetMethod + +composite type for the primal trust region method presented in [OberlinandWright-2006](@cite) +Contains the following fields: +-`silent`: bool setting the optimization solver used to silent (true) or not (false) +-`tol`: tolerance of the method +-`linear_program_solver`: linear solver used for the LP subproblem solved in the method +-`solver_attributes`: additional JuMP options to be given to the solver +-`ν`: penalty parameter, see [OberlinandWright-2006](@cite) +-`Δ`: trust region parameter, see [OberlinandWright-2006](@cite) + +# Reference +[OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming + +""" + +struct ActiveMethodLP_P <: AbstractActiveSetMethod + silent::Bool + tol::Float64 + linear_program_solver::DataType + solver_attributes::Tuple{Vararg{Pair}} + ν::Float64 + Δ::Float64 +end + +""" + ActiveMethodeMethodLP_P(linear_program_solver; kwargs...) + +Create an ActiveMethodeMethodLPEC where all fields can be specified as keyword arguments, a linear program solver must be provided. +The following default values are set +-`silent`: true +-`tol`: 1e-8 +-`solver_attributes`: () i.e. none +-`ν`: NaN which is replaced by 2*max(norm(λ, Inf), 1) when the associated find_active() is used, with λ the multipliers of the problem +-`Δ`:-NaN which is replaced by 4/(m+n) when the associated find_active() is used, with n the number of variables and m of constraints +""" + +ActiveMethodLP_P(linear_program_solver::DataType ;silent = true, tol = 1e-8, solver_attributes = (), ν = NaN, Δ = NaN) = ActiveMethodLP_P(silent, tol, linear_program_solver, solver_attributes, ν,Δ) + + +""" + ActiveMethodLPEC <: AbstractActiveSetMethod + +composite type for the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite) +Contains the following fields: +-`silent`: bool setting the optimization solver used to silent (true) or not (false) +-`tol`: tolerance of the method +-`solver`: solver used for the MILP subproblem solved in the method +-`solver_attributes`: additional JuMP options to be given to the solver +-`M`: big M constant, see [OberlinandWright-2006](@cite) +-`β`: positive test parameter, see [OberlinandWright-2006](@cite) +-`σ`: test parameter in (0,1), see [OberlinandWright-2006](@cite) + +# Reference +[OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming + +""" + +struct ActiveMethodLPEC <: AbstractActiveSetMethod + silent::Bool + tol::Float64 + solver::DataType + solver_attributes::Tuple{Vararg{Pair}} + M::Float64 + β::Float64 + σ::Float64 +end + +""" + ActiveMethodeMethodLPEC(solver; kwargs...) + +Create an ActiveMethodeMethodLPEC where all fields can be specified as keyword arguments, a solver must be provided. +The following default values are set: +-`silent`: true +-`tol`: 1e-8 +-`solver_attributes`: () i.e. none +-`M`: NaN wich is replaced by 10*max(norm(λ_ineq, Inf), norm(c_ineq, Inf)) when the associated find_active() is used, with c_ineq is the inequality constraints and λ_ineq their associated multipliers +-`β`: NaN wich is replaced by 1/(m+n) when the associated find_active() is used, with n the number of variables and m of constraints +-`σ`: 0.75 +""" + +ActiveMethodLPEC(solver::DataType ;silent = true, tol = 1e-8, solver_attributes = (), M = NaN, β = NaN, σ = 0.75) = ActiveMethodLPEC(silent, tol, solver, solver_attributes, M, β, σ) + +""" + ActiveMethodLPEC_A <: AbstractActiveSetMethod + +composite type for the linear programming approximation to the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite) +Contains the following fields: +-`silent`: bool setting the optimization solver used to silent (true) or not (false) +-`tol`: tolerance of the method +-`linear_program_solver`: linear solver used for the LP subproblem solved in the method +-`solver_attributes`: additional JuMP options to be given to the solver +-`β`: positive test parameter, see [OberlinandWright-2006](@cite) +-`σ`: test parameter in (0,1), see [OberlinandWright-2006](@cite) + +# Reference +[OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming + +""" + +struct ActiveMethodLPEC_A <: AbstractActiveSetMethod + silent::Bool + tol::Float64 + linear_program_solver::DataType + solver_attributes::Tuple{Vararg{Pair}} + β::Float64 + σ::Float64 +end + +""" +ActiveMethodeMethodLPEC_A(linear_program_solver; kwargs...) + +Create an ActiveMethodeMethodLPEC_A where all fields can be specified as keyword arguments, linear_program_solver must be provided. +The following default values are set: +-`silent`: true +-`tol`: 1e-8 +-`solver_attributes`: () i.e. none +-`β`: NaN wich is replaced by 1/(m+n) when the associated find_active() is used, with n the number of variables and m of constraints +-`σ`: 0.90 +""" + +ActiveMethodLPEC_A(linear_program_solver::DataType ;silent = true, tol = 1e-8, solver_attributes = (), β = NaN, σ = 0.90) = ActiveMethodLPEC_A(silent, tol, linear_program_solver, solver_attributes, β, σ) + + +""" + find_active(nlp, results, method::ActiveMethodeNaive) + +Finds the active set at results using the simple test c(x) = 0 + +Returns named tuple with 2 attributes: +- `active` : active set of usual constraints found +- `active_boundary` : active set of boundary constraints found + +""" + +function find_active(nlp, results, method::ActiveMethodSimple) + tol = method.tol + + x = results.solution + + constraints = NLPModels.cons(nlp, x) + + lvar = nlp.meta.lvar + uvar = nlp.meta.uvar + + lcon = nlp.meta.lcon + ucon = nlp.meta.ucon + + jlow = nlp.meta.jlow + jupp = nlp.meta.jupp + jrng = nlp.meta.jrng + + ilow = nlp.meta.ilow + iupp = nlp.meta.iupp + irng = nlp.meta.irng + + + var_ineq_idx = vcat(irng, iupp, ilow) + con_ineq_idx = vcat(jrng, jlow, jupp) + + active = Int64[] + active_boundary = Int64[] + + # Test run + for i in var_ineq_idx + if x[i] - lvar[i] <= tol || uvar[i] - x[i] <= tol + push!(active_boundary, i) + end + end + + for i in con_ineq_idx + if constraints[i] - lcon[i] <= tol || ucon[i] - constraints[i] <= tol + push!(active, i) + end + end + + return (active = active, active_boundary = active_boundary) +end + +""" + find_active(nlp, results, method::ActiveMethodeMethodLP_P) + +Implement the primal trust region method presented in [OberlinandWright-2006](@cite). For finding the active set at a solution with results being a point near the solution. + +Returns named tuple with 2 attributes: +- `active` : active set of usual constraints found +- `active_boundary` : active set of boundary constraints found + +# Reference +[OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming +""" + + +function find_active(nlp, results, method::ActiveMethodLP_P) + n = NLPModels.get_nvar(nlp) + m = NLPModels.get_ncon(nlp) + tol = method.tol + + jfix = nlp.meta.jfix + jlow = nlp.meta.jlow + jupp = nlp.meta.jupp + jrng = nlp.meta.jrng + + ifix = nlp.meta.ifix + ilow = nlp.meta.ilow + iupp = nlp.meta.iupp + irng = nlp.meta.irng + + allupp_con = vcat(jrng, jupp) + alllow_con = vcat(jrng, jlow) + + allupp_var = vcat(irng, iupp) + alllow_var = vcat(irng, ilow) + + lvar = nlp.meta.lvar + uvar = nlp.meta.uvar + + lcon = nlp.meta.lcon + ucon = nlp.meta.ucon + + x = results.solution + + y = results.multipliers + zl = results.multipliers_L + zu = results.multipliers_U + + Ji, Jj = NLPModels.jac_structure(nlp) + Jx = NLPModels.jac_coord(nlp, results.solution) + Jac = sparse(Ji, Jj, Jx, m, n) + + constraints = NLPModels.cons(nlp, x) + + g = NLPModels.grad(nlp, x) + + # Creating equality and inequality constraints and gradients in the right format + c_eq = vcat(constraints[jfix] - ucon[jfix], x[ifix] - uvar[ifix]) + n_eq = length(c_eq) + J_eq = vcat(Jac[jfix,:], spdiagm(ones(n))[ifix,:]) + + # we double the range constraints and use the convention c(x) <= 0 + c_ineq = vcat(constraints[allupp_con] - ucon[allupp_con], + lcon[alllow_con] - constraints[alllow_con], + x[allupp_var] - uvar[allupp_var], + lvar[alllow_var] - x[alllow_var] + ) + indices_to_constraints = vcat(allupp_con, alllow_con, allupp_var, alllow_var) + n_ineq = length(c_ineq) + J_ineq = vcat(Jac[allupp_con,:], -Jac[alllow_con,:], spdiagm(ones(n))[allupp_var,:], -spdiagm(ones(n))[alllow_var,:]) + + + # Parameter calculation + if isnan(method.ν) + ν = 2 * max(norm(y,Inf), norm(zl, Inf), norm(zu, Inf), 1) + else + ν = method.ν + end + if isnan(method.Δ) + Δ = 4/n + else + Δ = method.Δ + end + + # Error testing + if Δ <= 0 + throw(DomainError(Δ, "Δ must be positive")) + end + if ν <= 0 + throw(DomainError(ν, "ν must be positive")) + end + + # Sub problem resolution + model = Model(optimizer_with_attributes(method.linear_program_solver, method.solver_attributes...)) + + if method.silent + set_silent(model) + end + + @variable(model, -Δ <= d[1:n] <= Δ) # Trust region constraint + + if n_eq != 0 + @variable(model, u[1:n_eq] >= 0) + @variable(model, v[1:n_eq] >= 0) + + @constraint(model, J_eq * d + c_eq == u - v) + + @expression(model, slackCost_eq, ν * (sum(u)+sum(v))) + else + @expression(model, slackCost_eq, 0.) + end + + if n_ineq != 0 + @variable(model, r[1:n_ineq] >= 0) + + @constraint(model, J_ineq * d + c_ineq .<= r) + + @expression(model, slackCost_ineq, ν * sum(r)) + else + @expression(model, slackCost_eq, 0.) + end + + @objective(model, Min, dot(g,d) + slackCost_eq + slackCost_ineq) + + optimize!(model) + + if termination_status(model) != OPTIMAL + error("LP subproblem failed to converge : $(termination_status(model))") + end + + d_sol = value.(d) + + # Active set test + + active = Int64[] + active_boundary = Int64[] + + for i in 1:n_ineq + if dot(J_ineq[i, :], d_sol) + c_ineq[i] >= -tol + if i <= length(allupp_con) + length(alllow_con) + !(i in active) && push!(active, indices_to_constraints[i]) + else + !(i in active_boundary) && push!(active_boundary, indices_to_constraints[i]) + end + end + end + + + return (active = active, active_boundary = active_boundary) +end + +""" + find_active(nlp, results, method::ActiveMethodeMethodLPEC) + + +Implements the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite). For finding the active set at a solution with results being a point near the solution. + +Returns named tuple with 2 attributes: +- `active` : active set of usual constraints found +- `active_boundary` : active set of boundary constraints found + +# Reference +[OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming +""" + +function find_active(nlp, results, method::ActiveMethodLPEC) + n = NLPModels.get_nvar(nlp) + m = NLPModels.get_ncon(nlp) + tol = method.tol + + jfix = nlp.meta.jfix + jlow = nlp.meta.jlow + jupp = nlp.meta.jupp + jrng = nlp.meta.jrng + + ifix = nlp.meta.ifix + ilow = nlp.meta.ilow + iupp = nlp.meta.iupp + irng = nlp.meta.irng + + allupp_con = vcat(jrng, jupp) + alllow_con = vcat(jrng, jlow) + + allupp_var = vcat(irng, iupp) + alllow_var = vcat(irng, ilow) + + lvar = nlp.meta.lvar + uvar = nlp.meta.uvar + + lcon = nlp.meta.lcon + ucon = nlp.meta.ucon + + x = results.solution + + y = results.multipliers + zl = results.multipliers_L + zu = results.multipliers_U + + Ji, Jj = NLPModels.jac_structure(nlp) + Jx = NLPModels.jac_coord(nlp, x) + Jac = sparse(Ji, Jj, Jx, m, n) + + constraints = NLPModels.cons(nlp, x) + g = NLPModels.grad(nlp, x) + + # Creating equality and inequality constraints and gradients in the right format + c_eq = vcat(constraints[jfix] - ucon[jfix], x[ifix] - uvar[ifix]) + n_eq = length(c_eq) + J_eq = vcat(Jac[jfix,:], spdiagm(ones(n))[ifix,:]) + Jt_eq = transpose(J_eq) + + # we double the range constraints and use the convention c(x) <= 0 + c_ineq = vcat(constraints[allupp_con] - ucon[allupp_con], + lcon[alllow_con] - constraints[alllow_con], + x[allupp_var] - uvar[allupp_var], + lvar[alllow_var] - x[alllow_var] + ) + indices_to_constraints = vcat(allupp_con, alllow_con, allupp_var, alllow_var) + n_ineq = length(c_ineq) + J_ineq = vcat(Jac[allupp_con,:], -Jac[alllow_con,:], spdiagm(ones(n))[allupp_var,:], -spdiagm(ones(n))[alllow_var,:]) + Jt_ineq = transpose(J_ineq) + + # Parameter calculation + if isnan(method.M) + M = 10* max(norm(c_ineq, Inf), norm(y, Inf), norm(zu, Inf), norm(zl, Inf)) + else + M = method.M + end + if isnan(method.β) + β = 1/(m+n) + else + β = method.β + end + + σ = method.σ + + # Error testing + if M <= 0 + throw(DomainError(M, " M must be positive")) + end + if β <= 0 + throw(DomainError(β, " β must be positive")) + end + if σ <= 0 || σ >=1 + throw(DomainError(σ, "σ must be in (0,1)")) + end + + + # Sub problem resolution + model = Model(optimizer_with_attributes(method.solver, method.solver_attributes...)) + + if method.silent + set_silent(model) + end + + @variable(model, u[1:n]>=0) + @variable(model, v[1:n]>=0) + + if n_eq != 0 + @variable(model, λ_eq[1:n_eq]) + @expression(model, Jtprod_eq, Jt_eq * λ_eq) + else + @expression(model, Jtprod_eq, zeros(n)) + end + + if n_ineq != 0 + @variable(model, λ_ineq[1:n_ineq] >= 0) + @variable(model, s[1:n_ineq]) + @variable(model, y[1:n_ineq], Bin) + set_lower_bound.(s, max.(c_ineq, zeros(n_ineq))) + + @constraint(model, -c_ineq - s .<= -c_ineq .* y) + @constraint(model, λ_ineq - s .<= M*(1 .- y)) + + + @expression(model, Jtprod_ineq, Jt_ineq * λ_ineq) + @expression(model, slackCost, sum(s)) + else + @expression(model, Jtprod_ineq, zeros(n)) + @expression(model, slackCost, 0.) + end + + @constraint(model, Jtprod_eq + Jtprod_ineq + g == u - v) + + @objective(model, Min, slackCost + sum(u) + sum(v)) + + optimize!(model) + + if termination_status(model) != OPTIMAL + error("LP subproblem failed to converge : $(termination_status(model))") + end + + ω = objective_value(model) + norm(c_eq,1) + + if ω < -tol + trow(DomainError(ω, "Optimal value of subproblem is negative")) + end + + clamp(ω, 0., Inf) # TODO better way to deal with the case where ω is in [-tol, 0) ? + + # Active set test + + active = Int64[] + active_boundary = Int64[] + + for i in 1:n_ineq + if c_ineq[i] >= -(β*ω)^σ - tol + if i <= length(allupp_con) + length(alllow_con) + !(i in active) && push!(active, indices_to_constraints[i]) + else + !(i in active_boundary) && push!(active_boundary, indices_to_constraints[i]) + end + end + end + + return (active = active, active_boundary = active_boundary) +end + + +""" + find_active(nlp, results, method::ActiveMethodeMethodLPEC_A) + +Implements the linear programming approximation of the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite). For finding the active set at a solution with results being a point near the solution. + +Returns named tuple with 2 attributes: +- `active` : active set of usual constraints found +- `active_boundary` : active set of boundary constraints found + +# Reference +[OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming +""" + +function find_active(nlp, results, method::ActiveMethodLPEC_A) + n = NLPModels.get_nvar(nlp) + m = NLPModels.get_ncon(nlp) + tol = method.tol + + jfix = nlp.meta.jfix + jlow = nlp.meta.jlow + jupp = nlp.meta.jupp + jrng = nlp.meta.jrng + + ifix = nlp.meta.ifix + ilow = nlp.meta.ilow + iupp = nlp.meta.iupp + irng = nlp.meta.irng + + allupp_con = vcat(jrng, jupp) + alllow_con = vcat(jrng, jlow) + + allupp_var = vcat(irng, iupp) + alllow_var = vcat(irng, ilow) + + lvar = nlp.meta.lvar + uvar = nlp.meta.uvar + + lcon = nlp.meta.lcon + ucon = nlp.meta.ucon + + + x = results.solution + y = results.multipliers + zl = results.multipliers_L + zu = results.multipliers_U + + Ji, Jj = NLPModels.jac_structure(nlp) + Jx = NLPModels.jac_coord(nlp, results.solution) + Jac = sparse(Ji, Jj, Jx, m, n) + constraints = NLPModels.cons(nlp, x) + + g = NLPModels.grad(nlp, x) + + # Creating equality and inequality constraints and gradients in the right format + c_eq = vcat(constraints[jfix] - ucon[jfix], x[ifix] - uvar[ifix]) + n_eq = length(c_eq) + J_eq = vcat(Jac[jfix,:], spdiagm(ones(n))[ifix,:]) + Jt_eq = transpose(J_eq) + + # we double the range constraints and use the convention c(x) <= 0 + c_ineq = vcat(constraints[allupp_con] - ucon[allupp_con], + lcon[alllow_con] - constraints[alllow_con], + x[allupp_var] - uvar[allupp_var], + lvar[alllow_var] - x[alllow_var] + ) + indices_to_constraints = vcat(allupp_con, alllow_con, allupp_var, alllow_var) + n_ineq = length(c_ineq) + neg_idx = findall(x -> x<-tol, c_ineq) + pos_idx = setdiff(1:n_ineq, neg_idx) + J_ineq = vcat(Jac[allupp_con,:], -Jac[alllow_con,:], spdiagm(ones(n))[allupp_var,:], -spdiagm(ones(n))[alllow_var,:]) + Jt_ineq = transpose(J_ineq) + + # Parameter calculation + K_1 = max(norm(c_ineq, Inf), norm(y, Inf), norm(zl, Inf), norm(zu, Inf)) + + if isnan(method.β) + β = 1/(m+n) + else + β = method.β + end + σ = method.σ + + # Error testing + if β <= 0 + throw(DomainError(β, " β must be positive")) + end + if σ <= 0 || σ >=1 + throw(DomainError(σ, "σ must be in (0,1)")) + end + + # Sub problem resolution + model = Model(optimizer_with_attributes(method.linear_program_solver, method.solver_attributes...)) + + if method.silent + set_silent(model) + end + + + @variable(model, u[1:n]>=0) + @variable(model, v[1:n]>=0) + + if n_eq != 0 + @variable(model, λ_eq[1:n_eq]) + @expression(model, Jtprod_eq, Jt_eq * λ_eq) + else + @expression(model, Jtprod_eq, zeros(n)) + end + + if n_ineq != 0 + @variable(model, K_1 >= λ_ineq[1:n_ineq] >= 0) + + @expression(model, Jtprod_ineq, Jt_ineq * λ_ineq) + @expression(model, ineqCost, -dot(c_ineq[neg_idx], λ_ineq[neg_idx])) + else + @expression(model, Jtprod_ineq, zeros(n)) + @expression(model, ineqCost, 0.) + end + + @constraint(model, Jtprod_eq + Jtprod_ineq + g == u - v) + + @objective(model, Min, ineqCost + sum(u) + sum(v)) + + optimize!(model) + + if termination_status(model) != OPTIMAL + error("LP subproblem failed to converge : $(termination_status(model))") + end + + if n_eq != 0 + λ_eq_sol = value.(λ_eq) + else + λ_eq_sol = Float64[] + end + + if n_ineq != 0 + λ_ineq_sol = value.(λ_ineq) + else + λ_ineq_sol = Float64[] + end + + #TODO find a better solution than clamp() to ensure ρ_sup is positive (under tol) + ρ_sup = sum(clamp(-c_ineq[i] * λ_ineq_sol[i], 0., Inf)^(1/2) for i in neg_idx; init=0.) + + sum(c_ineq[pos_idx]) + + norm(c_eq, 1) + + norm(Jt_ineq * λ_ineq_sol + Jt_eq * λ_eq_sol + g, 1) + + + if ρ_sup < -tol + trow(DomainError(ρ_sup, "test bound ρ_sup is negative")) + end + ρ_sup = clamp(ρ_sup, 0., Inf) #TODO same question here + + # Active set test + + active = Int64[] + active_boundary = Int64[] + for i in 1:n_ineq + if c_ineq[i] >= -(β*ρ_sup)^σ - tol + if i <= length(allupp_con) + length(alllow_con) + !(i in active) && push!(active, indices_to_constraints[i]) + else + !(i in active_boundary) && push!(active_boundary, indices_to_constraints[i]) + end + end + end + + return (active = active, active_boundary = active_boundary) +end \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml index 17dc64c..b47710f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,3 +5,4 @@ MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsJuMP = "792afdf1-32c1-5681-94e0-d7bf7a5df49e" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index dac0ace..1d3310c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test using MadNLP using NLPModels using MadCheck +using HiGHS include("instances.jl") @@ -17,3 +18,26 @@ include("instances.jl") @test chk.inf_cc_con <= 1e-8 end +@testset "Test active set detection" begin + nlp = hs15_model() + results = madnlp(nlp; print_level=MadNLP.ERROR) + + n = NLPModels.get_nvar(nlp) + + new_point = results + max_dist = 0 + for i in 1:n + noise = 2*rand() - 1 + new_point.solution[i] += (max_dist*noise)/n + end + + active_solution = MadCheck.find_active(nlp, results,MadCheck.ActiveMethodSimple()) + + active_LP_P = MadCheck.find_active(nlp, new_point, MadCheck.ActiveMethodLP_P(HiGHS.Optimizer; Δ = 4 * 1e-8/n)) + active_LPEC = MadCheck.find_active(nlp, new_point, MadCheck.ActiveMethodLPEC(HiGHS.Optimizer)) + active_LPEC_A = MadCheck.find_active(nlp, new_point, MadCheck.ActiveMethodLPEC_A(HiGHS.Optimizer)) + + @test Set(active_LP_P.active) == Set(active_solution.active) && Set(active_LP_P.active_boundary) == Set(active_solution.active_boundary) + @test Set(active_LPEC.active) == Set(active_solution.active) && Set(active_LPEC.active_boundary) == Set(active_solution.active_boundary) + @test Set(active_LPEC_A.active) == Set(active_solution.active) && Set(active_LPEC_A.active_boundary) == Set(active_solution.active_boundary) +end \ No newline at end of file From 0b66659c7a234d3ef698ce58466ddbb781a8921e Mon Sep 17 00:00:00 2001 From: GhislainDv Date: Thu, 27 Aug 2026 15:59:08 +0200 Subject: [PATCH 2/2] Made changes to active set & Implemented LICQ and the SVD method --- src/MadCheck.jl | 3 +- src/active_set.jl | 251 +++++++++++++++++++------------- src/constraint_qualification.jl | 91 +++++++++++- src/jacobian_degeneracy.jl | 66 +++++++++ test/instances.jl | 29 ++++ test/runtests.jl | 23 ++- 6 files changed, 351 insertions(+), 112 deletions(-) create mode 100644 src/jacobian_degeneracy.jl diff --git a/src/MadCheck.jl b/src/MadCheck.jl index ebe3c7e..dd00944 100644 --- a/src/MadCheck.jl +++ b/src/MadCheck.jl @@ -6,8 +6,9 @@ import NLPModels using JuMP include("solution.jl") -include("constraint_qualification.jl") include("active_set.jl") +include("jacobian_degeneracy.jl") +include("constraint_qualification.jl") include("feasibility.jl") end # module MadCheck diff --git a/src/active_set.jl b/src/active_set.jl index 9f477e6..c77bfcf 100644 --- a/src/active_set.jl +++ b/src/active_set.jl @@ -3,39 +3,36 @@ Base type for the different active set methods implemented """ - abstract type AbstractActiveSetMethod end """ - ActiveMethodSimple <: AbstractActiveSetMethod + BasicActiveSet <: AbstractActiveSetMethod -composite type for the naive method consisting of runing the test c(x) = 0, +Composite type for the basic method consisting of runing the test c(x) = 0, contains one parameter `tol`: tolerance of the method """ - -struct ActiveMethodSimple <: AbstractActiveSetMethod +struct BasicActiveSet <: AbstractActiveSetMethod tol::Float64 end """ - ActiveMethodSimple(; kwargs...) + BasicActiveSet(; kwargs...) Creates an ActiveMethodeSimple with `tol` = 1e-8. """ - -ActiveMethodSimple(;tol = 1e-8) = ActiveMethodSimple(tol) +BasicActiveSet(;tol = 1e-8) = BasicActiveSet(tol) """ - ActiveMethodLP_P <: AbstractActiveSetMethod + PrimalActiveSetLP <: AbstractActiveSetMethod -composite type for the primal trust region method presented in [OberlinandWright-2006](@cite) +Composite type for the primal trust region method presented in [OberlinandWright-2006](@cite). We note that this method cannot detect weakly active constraints. Contains the following fields: -`silent`: bool setting the optimization solver used to silent (true) or not (false) -`tol`: tolerance of the method -`linear_program_solver`: linear solver used for the LP subproblem solved in the method --`solver_attributes`: additional JuMP options to be given to the solver +-`solver_options`: additional JuMP options to be given to the solver -`ν`: penalty parameter, see [OberlinandWright-2006](@cite) -`Δ`: trust region parameter, see [OberlinandWright-2006](@cite) @@ -43,40 +40,37 @@ Contains the following fields: [OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming """ - -struct ActiveMethodLP_P <: AbstractActiveSetMethod +struct PrimalActiveSetLP <: AbstractActiveSetMethod silent::Bool tol::Float64 linear_program_solver::DataType - solver_attributes::Tuple{Vararg{Pair}} + solver_options::Tuple{Vararg{Pair}} ν::Float64 Δ::Float64 end """ - ActiveMethodeMethodLP_P(linear_program_solver; kwargs...) + PrimalActiveSetLP(linear_program_solver; kwargs...) -Create an ActiveMethodeMethodLPEC where all fields can be specified as keyword arguments, a linear program solver must be provided. +Create an PrimalActiveSetLP where all fields can be specified as keyword arguments, a linear program solver and the trust region parameter must be provided. The following default values are set -`silent`: true -`tol`: 1e-8 --`solver_attributes`: () i.e. none +-`solver_options`: () i.e. none -`ν`: NaN which is replaced by 2*max(norm(λ, Inf), 1) when the associated find_active() is used, with λ the multipliers of the problem --`Δ`:-NaN which is replaced by 4/(m+n) when the associated find_active() is used, with n the number of variables and m of constraints """ - -ActiveMethodLP_P(linear_program_solver::DataType ;silent = true, tol = 1e-8, solver_attributes = (), ν = NaN, Δ = NaN) = ActiveMethodLP_P(silent, tol, linear_program_solver, solver_attributes, ν,Δ) +PrimalActiveSetLP(linear_program_solver::DataType, Δ;silent = true, tol = 1e-8, solver_options = (), ν = NaN) = PrimalActiveSetLP(silent, tol, linear_program_solver, solver_options, ν,Δ) """ - ActiveMethodLPEC <: AbstractActiveSetMethod + PrimalDualActiveSetLPEC <: AbstractActiveSetMethod -composite type for the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite) +Composite type for the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite) Contains the following fields: -`silent`: bool setting the optimization solver used to silent (true) or not (false) -`tol`: tolerance of the method -`solver`: solver used for the MILP subproblem solved in the method --`solver_attributes`: additional JuMP options to be given to the solver +-`solver_options`: additional JuMP options to be given to the solver -`M`: big M constant, see [OberlinandWright-2006](@cite) -`β`: positive test parameter, see [OberlinandWright-2006](@cite) -`σ`: test parameter in (0,1), see [OberlinandWright-2006](@cite) @@ -86,40 +80,40 @@ Contains the following fields: """ -struct ActiveMethodLPEC <: AbstractActiveSetMethod +struct PrimalDualActiveSetLPEC <: AbstractActiveSetMethod silent::Bool tol::Float64 solver::DataType - solver_attributes::Tuple{Vararg{Pair}} + solver_options::Tuple{Vararg{Pair}} M::Float64 β::Float64 σ::Float64 end """ - ActiveMethodeMethodLPEC(solver; kwargs...) + PrimalDualActiveSetLPEC(solver; kwargs...) -Create an ActiveMethodeMethodLPEC where all fields can be specified as keyword arguments, a solver must be provided. +Create an PrimalDualActiveSetLPEC where all fields can be specified as keyword arguments, a MILP optimization solver must be provided. The following default values are set: -`silent`: true -`tol`: 1e-8 --`solver_attributes`: () i.e. none +-`solver_options`: () i.e. none -`M`: NaN wich is replaced by 10*max(norm(λ_ineq, Inf), norm(c_ineq, Inf)) when the associated find_active() is used, with c_ineq is the inequality constraints and λ_ineq their associated multipliers -`β`: NaN wich is replaced by 1/(m+n) when the associated find_active() is used, with n the number of variables and m of constraints -`σ`: 0.75 """ -ActiveMethodLPEC(solver::DataType ;silent = true, tol = 1e-8, solver_attributes = (), M = NaN, β = NaN, σ = 0.75) = ActiveMethodLPEC(silent, tol, solver, solver_attributes, M, β, σ) +PrimalDualActiveSetLPEC(solver::DataType ;silent = true, tol = 1e-8, solver_options = (), M = NaN, β = NaN, σ = 0.75) = PrimalDualActiveSetLPEC(silent, tol, solver, solver_options, M, β, σ) """ - ActiveMethodLPEC_A <: AbstractActiveSetMethod + ApproximatePrimalDualActiveSetLPEC <: AbstractActiveSetMethod -composite type for the linear programming approximation to the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite) +Composite type for the linear programming approximation to the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite) Contains the following fields: -`silent`: bool setting the optimization solver used to silent (true) or not (false) -`tol`: tolerance of the method -`linear_program_solver`: linear solver used for the LP subproblem solved in the method --`solver_attributes`: additional JuMP options to be given to the solver +-`solver_options`: additional JuMP options to be given to the solver -`β`: positive test parameter, see [OberlinandWright-2006](@cite) -`σ`: test parameter in (0,1), see [OberlinandWright-2006](@cite) @@ -128,42 +122,41 @@ Contains the following fields: """ -struct ActiveMethodLPEC_A <: AbstractActiveSetMethod +struct ApproximatePrimalDualActiveSetLPEC <: AbstractActiveSetMethod silent::Bool tol::Float64 linear_program_solver::DataType - solver_attributes::Tuple{Vararg{Pair}} + solver_options::Tuple{Vararg{Pair}} β::Float64 σ::Float64 end """ -ActiveMethodeMethodLPEC_A(linear_program_solver; kwargs...) +ApproximatePrimalDualActiveSetLPEC(linear_program_solver; kwargs...) -Create an ActiveMethodeMethodLPEC_A where all fields can be specified as keyword arguments, linear_program_solver must be provided. +Create an ApproximatePrimalDualActiveSetLPEC where all fields can be specified as keyword arguments, linear_program_solver must be provided. The following default values are set: -`silent`: true -`tol`: 1e-8 --`solver_attributes`: () i.e. none +-`solver_options`: () i.e. none -`β`: NaN wich is replaced by 1/(m+n) when the associated find_active() is used, with n the number of variables and m of constraints -`σ`: 0.90 """ -ActiveMethodLPEC_A(linear_program_solver::DataType ;silent = true, tol = 1e-8, solver_attributes = (), β = NaN, σ = 0.90) = ActiveMethodLPEC_A(silent, tol, linear_program_solver, solver_attributes, β, σ) +ApproximatePrimalDualActiveSetLPEC(linear_program_solver::DataType ;silent = true, tol = 1e-8, solver_options = (), β = NaN, σ = 0.90) = ApproximatePrimalDualActiveSetLPEC(silent, tol, linear_program_solver, solver_options, β, σ) """ - find_active(nlp, results, method::ActiveMethodeNaive) - -Finds the active set at results using the simple test c(x) = 0 + find_active(nlp, results, method::BasicActiveSet) +0 +Find the active inequality constraints ``c_i`` matching their lower-bound (``c_i(x) = lb_i``) or their upper-bound (``c_i(x) = ub_i``) at the current solution ``x`` store in `results.solution` Returns named tuple with 2 attributes: - `active` : active set of usual constraints found - `active_boundary` : active set of boundary constraints found """ - -function find_active(nlp, results, method::ActiveMethodSimple) +function find_active(nlp, results, method::BasicActiveSet) tol = method.tol x = results.solution @@ -208,9 +201,10 @@ function find_active(nlp, results, method::ActiveMethodSimple) end """ - find_active(nlp, results, method::ActiveMethodeMethodLP_P) + find_active(nlp, results, method::PrimalActiveSetLP) Implement the primal trust region method presented in [OberlinandWright-2006](@cite). For finding the active set at a solution with results being a point near the solution. +We note that this method cannot detect weakly active constraints. Returns named tuple with 2 attributes: - `active` : active set of usual constraints found @@ -219,9 +213,7 @@ Returns named tuple with 2 attributes: # Reference [OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming """ - - -function find_active(nlp, results, method::ActiveMethodLP_P) +function find_active(nlp, results, method::PrimalActiveSetLP) n = NLPModels.get_nvar(nlp) m = NLPModels.get_ncon(nlp) tol = method.tol @@ -236,12 +228,6 @@ function find_active(nlp, results, method::ActiveMethodLP_P) iupp = nlp.meta.iupp irng = nlp.meta.irng - allupp_con = vcat(jrng, jupp) - alllow_con = vcat(jrng, jlow) - - allupp_var = vcat(irng, iupp) - alllow_var = vcat(irng, ilow) - lvar = nlp.meta.lvar uvar = nlp.meta.uvar @@ -267,16 +253,41 @@ function find_active(nlp, results, method::ActiveMethodLP_P) n_eq = length(c_eq) J_eq = vcat(Jac[jfix,:], spdiagm(ones(n))[ifix,:]) - # we double the range constraints and use the convention c(x) <= 0 - c_ineq = vcat(constraints[allupp_con] - ucon[allupp_con], - lcon[alllow_con] - constraints[alllow_con], - x[allupp_var] - uvar[allupp_var], - lvar[alllow_var] - x[alllow_var] + c_ineq = vcat(constraints[jrng], + constraints[jupp], + constraints[jlow], + x[irng], + x[iupp], + x[ilow], + ) + + u = vcat(ucon[jrng], + ucon[jupp], + ucon[jlow], + uvar[irng], + uvar[iupp], + uvar[ilow], ) - indices_to_constraints = vcat(allupp_con, alllow_con, allupp_var, alllow_var) - n_ineq = length(c_ineq) - J_ineq = vcat(Jac[allupp_con,:], -Jac[alllow_con,:], spdiagm(ones(n))[allupp_var,:], -spdiagm(ones(n))[alllow_var,:]) + l = vcat(lcon[jrng], + lcon[jupp], + lcon[jlow], + lvar[irng], + lvar[iupp], + lvar[ilow], + ) + + J_ineq = vcat(Jac[jrng, :], + Jac[jupp, :], + Jac[jlow, :], + spdiagm(ones(n))[irng,:], + spdiagm(ones(n))[iupp,:], + spdiagm(ones(n))[ilow,:], + ) + + indices_to_constraints = vcat(jrng, jupp, jlow, irng, iupp, ilow) + n_ineq = length(c_ineq) + n_ineq_con = length(jrng) + length(jupp) + length(jlow) # Parameter calculation if isnan(method.ν) @@ -284,11 +295,7 @@ function find_active(nlp, results, method::ActiveMethodLP_P) else ν = method.ν end - if isnan(method.Δ) - Δ = 4/n - else - Δ = method.Δ - end + Δ = method.Δ # Error testing if Δ <= 0 @@ -299,7 +306,7 @@ function find_active(nlp, results, method::ActiveMethodLP_P) end # Sub problem resolution - model = Model(optimizer_with_attributes(method.linear_program_solver, method.solver_attributes...)) + model = Model(optimizer_with_attributes(method.linear_program_solver, method.solver_options...)) if method.silent set_silent(model) @@ -320,8 +327,15 @@ function find_active(nlp, results, method::ActiveMethodLP_P) if n_ineq != 0 @variable(model, r[1:n_ineq] >= 0) - - @constraint(model, J_ineq * d + c_ineq .<= r) + + for i in 1:n_ineq + if isfinite(u[i]) + @constraint(model, J_ineq[i,:]' * d + c_ineq[i] <= u[i] + r[i]) + end + if isfinite(l[i]) + @constraint(model, J_ineq[i,:]' * d + c_ineq[i] >= l[i] + r[i]) + end + end @expression(model, slackCost_ineq, ν * sum(r)) else @@ -344,11 +358,11 @@ function find_active(nlp, results, method::ActiveMethodLP_P) active_boundary = Int64[] for i in 1:n_ineq - if dot(J_ineq[i, :], d_sol) + c_ineq[i] >= -tol - if i <= length(allupp_con) + length(alllow_con) - !(i in active) && push!(active, indices_to_constraints[i]) + if dot(J_ineq[i, :], d_sol) + c_ineq[i] >= u[i] - tol || dot(J_ineq[i, :], d_sol) + c_ineq[i] <= l[i] + tol + if i <= n_ineq_con + push!(active, indices_to_constraints[i]) else - !(i in active_boundary) && push!(active_boundary, indices_to_constraints[i]) + push!(active_boundary, indices_to_constraints[i]) end end end @@ -358,7 +372,7 @@ function find_active(nlp, results, method::ActiveMethodLP_P) end """ - find_active(nlp, results, method::ActiveMethodeMethodLPEC) + find_active(nlp, results, method::PrimalDualActiveSetLPEC) Implements the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite). For finding the active set at a solution with results being a point near the solution. @@ -370,8 +384,7 @@ Returns named tuple with 2 attributes: # Reference [OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming """ - -function find_active(nlp, results, method::ActiveMethodLPEC) +function find_active(nlp, results, method::PrimalDualActiveSetLPEC) n = NLPModels.get_nvar(nlp) m = NLPModels.get_ncon(nlp) tol = method.tol @@ -417,15 +430,30 @@ function find_active(nlp, results, method::ActiveMethodLPEC) J_eq = vcat(Jac[jfix,:], spdiagm(ones(n))[ifix,:]) Jt_eq = transpose(J_eq) - # we double the range constraints and use the convention c(x) <= 0 - c_ineq = vcat(constraints[allupp_con] - ucon[allupp_con], - lcon[alllow_con] - constraints[alllow_con], - x[allupp_var] - uvar[allupp_var], - lvar[alllow_var] - x[alllow_var] + # For inequality contraints, we use the convention c(x) <= 0 + + c_ineq = vcat(min.(constraints[jrng] - ucon[jrng], lcon[jrng] - constraints[jrng]), + constraints[jupp] - ucon[jupp], + lcon[jlow] - constraints[jlow], + min.(x[irng] - uvar[irng], lvar[irng] - x[irng]), + x[iupp] - uvar[iupp], + lvar[ilow] - x[ilow], + ) + + sign_con = [(ucon[i] - constraints[i]) > (constraints[i] - lcon[i]) ? -1 : 1 for i in jrng] + sign_var = [(uvar[i] - x[i]) > (x[i] - lvar[i]) ? -1 : 1 for i in jrng] + + J_ineq = vcat(Jac[jrng, :] .* sign_con, + Jac[jupp, :], + -Jac[jlow, :], + spdiagm(ones(n))[irng,:] .* sign_var, + spdiagm(ones(n))[iupp,:], + -spdiagm(ones(n))[ilow,:], ) - indices_to_constraints = vcat(allupp_con, alllow_con, allupp_var, alllow_var) + + indices_to_constraints = vcat(jrng, jupp, jlow, irng, iupp, ilow) n_ineq = length(c_ineq) - J_ineq = vcat(Jac[allupp_con,:], -Jac[alllow_con,:], spdiagm(ones(n))[allupp_var,:], -spdiagm(ones(n))[alllow_var,:]) + n_ineq_con = length(jrng) + length(jupp) + length(jlow) Jt_ineq = transpose(J_ineq) # Parameter calculation @@ -455,7 +483,7 @@ function find_active(nlp, results, method::ActiveMethodLPEC) # Sub problem resolution - model = Model(optimizer_with_attributes(method.solver, method.solver_attributes...)) + model = Model(optimizer_with_attributes(method.solver, method.solver_options...)) if method.silent set_silent(model) @@ -513,10 +541,10 @@ function find_active(nlp, results, method::ActiveMethodLPEC) for i in 1:n_ineq if c_ineq[i] >= -(β*ω)^σ - tol - if i <= length(allupp_con) + length(alllow_con) - !(i in active) && push!(active, indices_to_constraints[i]) + if i <= n_ineq_con + push!(active, indices_to_constraints[i]) else - !(i in active_boundary) && push!(active_boundary, indices_to_constraints[i]) + push!(active_boundary, indices_to_constraints[i]) end end end @@ -526,7 +554,7 @@ end """ - find_active(nlp, results, method::ActiveMethodeMethodLPEC_A) + find_active(nlp, results, method::ApproximatePrimalDualActiveSetLPEC) Implements the linear programming approximation of the method based on the primal-dual estimate presented in [OberlinandWright-2006](@cite). For finding the active set at a solution with results being a point near the solution. @@ -537,8 +565,7 @@ Returns named tuple with 2 attributes: # Reference [OberlinandWright-2006] Oberlin and Wright - 2006 - Active Set Identification in Nonlinear Programming """ - -function find_active(nlp, results, method::ActiveMethodLPEC_A) +function find_active(nlp, results, method::ApproximatePrimalDualActiveSetLPEC) n = NLPModels.get_nvar(nlp) m = NLPModels.get_ncon(nlp) tol = method.tol @@ -584,18 +611,34 @@ function find_active(nlp, results, method::ActiveMethodLPEC_A) J_eq = vcat(Jac[jfix,:], spdiagm(ones(n))[ifix,:]) Jt_eq = transpose(J_eq) - # we double the range constraints and use the convention c(x) <= 0 - c_ineq = vcat(constraints[allupp_con] - ucon[allupp_con], - lcon[alllow_con] - constraints[alllow_con], - x[allupp_var] - uvar[allupp_var], - lvar[alllow_var] - x[alllow_var] + # For inequality contraints, we use the convention c(x) <= 0 + + c_ineq = vcat(min.(constraints[jrng] - ucon[jrng], lcon[jrng] - constraints[jrng]), + constraints[jupp] - ucon[jupp], + lcon[jlow] - constraints[jlow], + min.(x[irng] - uvar[irng], lvar[irng] - x[irng]), + x[iupp] - uvar[iupp], + lvar[ilow] - x[ilow], + ) + + sign_con = [(ucon[i] - constraints[i]) > (constraints[i] - lcon[i]) ? -1 : 1 for i in jrng] + sign_var = [(uvar[i] - x[i]) > (x[i] - lvar[i]) ? -1 : 1 for i in jrng] + + J_ineq = vcat(Jac[jrng, :] .* sign_con, + Jac[jupp, :], + -Jac[jlow, :], + spdiagm(ones(n))[irng,:] .* sign_var, + spdiagm(ones(n))[iupp,:], + -spdiagm(ones(n))[ilow,:], ) - indices_to_constraints = vcat(allupp_con, alllow_con, allupp_var, alllow_var) + + indices_to_constraints = vcat(jrng, jupp, jlow, irng, iupp, ilow) n_ineq = length(c_ineq) + n_ineq_con = length(jrng) + length(jupp) + length(jlow) + Jt_ineq = transpose(J_ineq) + neg_idx = findall(x -> x<-tol, c_ineq) pos_idx = setdiff(1:n_ineq, neg_idx) - J_ineq = vcat(Jac[allupp_con,:], -Jac[alllow_con,:], spdiagm(ones(n))[allupp_var,:], -spdiagm(ones(n))[alllow_var,:]) - Jt_ineq = transpose(J_ineq) # Parameter calculation K_1 = max(norm(c_ineq, Inf), norm(y, Inf), norm(zl, Inf), norm(zu, Inf)) @@ -616,7 +659,7 @@ function find_active(nlp, results, method::ActiveMethodLPEC_A) end # Sub problem resolution - model = Model(optimizer_with_attributes(method.linear_program_solver, method.solver_attributes...)) + model = Model(optimizer_with_attributes(method.linear_program_solver, method.solver_options...)) if method.silent set_silent(model) @@ -669,8 +712,8 @@ function find_active(nlp, results, method::ActiveMethodLPEC_A) ρ_sup = sum(clamp(-c_ineq[i] * λ_ineq_sol[i], 0., Inf)^(1/2) for i in neg_idx; init=0.) + sum(c_ineq[pos_idx]) + norm(c_eq, 1) + - norm(Jt_ineq * λ_ineq_sol + Jt_eq * λ_eq_sol + g, 1) - + norm(Jt_ineq * λ_ineq_sol + Jt_eq * λ_eq_sol + g, 1 + ) if ρ_sup < -tol trow(DomainError(ρ_sup, "test bound ρ_sup is negative")) @@ -683,10 +726,10 @@ function find_active(nlp, results, method::ActiveMethodLPEC_A) active_boundary = Int64[] for i in 1:n_ineq if c_ineq[i] >= -(β*ρ_sup)^σ - tol - if i <= length(allupp_con) + length(alllow_con) - !(i in active) && push!(active, indices_to_constraints[i]) + if i <= n_ineq_con + push!(active, indices_to_constraints[i]) else - !(i in active_boundary) && push!(active_boundary, indices_to_constraints[i]) + push!(active_boundary, indices_to_constraints[i]) end end end diff --git a/src/constraint_qualification.jl b/src/constraint_qualification.jl index 0504da8..ddc9f1e 100644 --- a/src/constraint_qualification.jl +++ b/src/constraint_qualification.jl @@ -1,3 +1,90 @@ +# TODO add MFCQ -# TODO: -# check if LICQ, SCS or MFCQ are satisfied at current point +""" + build_work_jacobian(nlp, results, active, active_boundary) + +Create the jacobian the method is working on based on what active set wants to be used (i.e. active = [] and active_boundary = [] for MFCQ). This jacobian can then be given to methods in jacobian_degeneracy for analysis. +We set the sign convention c(x) <= 0 for use in MFCQ. + +Return +-`J_work`: Constructed jacobian +-`indices_to_constraints`: vector mapping the indices of the columns of `J_work` to the constraints indices, the first part of it maps to usual constraints, the second to bound contraints. +""" +function build_work_jacobian(nlp, results, active, active_boundary) + n = NLPModels.get_nvar(nlp) + m = NLPModels.get_ncon(nlp) + + jfix = nlp.meta.jfix + ifix = nlp.meta.ifix + + ucon = nlp.meta.ucon + lcon = nlp.meta.lcon + uvar = nlp.meta.lvar + lvar = nlp.meta.uvar + + x = results.solution + Ji, Jj = NLPModels.jac_structure(nlp) + Jx = NLPModels.jac_coord(nlp, x) + J = sparse(Ji, Jj, Jx, m, n) + + constraints = NLPModels.cons(nlp, x) + + sign = [(ucon[i] - constraints[i]) > (constraints[i] - lcon[i]) ? -1 : 1 for i in active] + sign_boundary = [(uvar[i] - x[i]) > (x[i] - lvar[i]) ? -1 : 1 for i in active_boundary] + + J_work = vcat( + J[jfix, :], + J[active, :] .* sign, + spdiagm(ones(n))[ifix, :], + spdiagm(ones(n))[active_boundary, :] .* sign_boundary, + ) + + indices_to_constraints = vcat(jfix, active, ifix, active_boundary) + + return J_work, indices_to_constraints +end + +""" + test_LICQ(nlp, results, active_method::AbstractActiveSetMethod, degen_method::AbstractDegenJacMethod) + +Tests the LICQ condition at a given point by looking for dependent constraits in the set of active and equality constraints. +Takes the following arguments: +-`nlp`: non linear programming problem studied +-`results`: Point studied +-`active_method`: method used for finding the active set +-`degen_method`: method used for finding degenerate constraints in the set of equality and active constraints. + +Returns a list of named tuples with 2 fields: `constraints` and `bounds` corresponding to a set of dependent constraints/bounds. +""" +function test_LICQ( + nlp, + results, + active_method::AbstractActiveSetMethod, + degen_method::AbstractDegenJacMethod, +) + + active, active_boundary = find_active(nlp, results, active_method) + + Jac, indices_to_constraints = build_work_jacobian(nlp, results, active, active_boundary) + + list_degen_cons = find_degenerate(Jac, degen_method) + + # Return + rep = [] + n_work_con = length(nlp.meta.jfix) + length(active) + for degen_cons in list_degen_cons + cons = [] + bounds = [] + for i in degen_cons + if i <= n_work_con + push!(cons, indices_to_constraints[i]) + else + push!(bounds, indices_to_constraints[i]) + end + end + + push!(rep, (constraints = cons, bounds = bounds)) + end + + return rep +end diff --git a/src/jacobian_degeneracy.jl b/src/jacobian_degeneracy.jl new file mode 100644 index 0000000..02606bc --- /dev/null +++ b/src/jacobian_degeneracy.jl @@ -0,0 +1,66 @@ +# TODO implement Degen Hunter, Qr and maybe dulmage-mendelsohn + +""" + AbstractDegenJacMethod + +Base type for the different methods implemented for determining sets of degenerate columns in a given Jacobian +""" +abstract type AbstractDegenJacMethod end + +""" + DegenJacSVD <: AbstractDegenJacMethod + +Composite type for the method based on the singular value decomposition, +contains one parameter + `tol`: tolerance of the method +""" +struct DegenJacSVD <: AbstractDegenJacMethod + tol::Float64 +end + + +""" + DegenJacSVD(; kwargs...) + +Creates an DegenJacSVD with `tol` = 1e-8. + +""" +DegenJacSVD(; tol = 1e-8) = DegenJacSVD(tol) + +""" + find_degenerate(Jac, method::DegenJacSVD) + +Finds the dependent columns in the matrice Jac using it's singular value decomposition. + +Return +list of vectors of indices corresponding to dependent sets of columns of Jac +""" +function find_degenerate(Jac, method::DegenJacSVD) + n_jac, n = size(Jac) + tol = method.tol + + if n_jac > n # If n_cons_eq > n we need to use the full svd to have access to all the columns of U to then run the test on all of them + J_svd = svd(Array(Jac), full = true) + else + J_svd = svd(Array(Jac)) + end + U = J_svd.U + S = J_svd.S + + r = min(n_jac, n) + # Rank + while r >= 1 && S[r] <= tol + r-=1 + end + + list_degen_colums = [] + + for u in eachcol(U[:, (r+1):n_jac]) + degen_columns = findall(x -> abs(x)>tol, u) + if !(Set(degen_columns) in Set.(list_degen_colums)) + push!(list_degen_colums, degen_columns) + end + end + + return list_degen_colums +end diff --git a/test/instances.jl b/test/instances.jl index 4bf8768..9e3d4a7 100644 --- a/test/instances.jl +++ b/test/instances.jl @@ -12,3 +12,32 @@ function hs15_model() @constraint(model, x[1] + x[2]^2 ≥ 0) return MathOptNLPModel(model) end + +function degenerate_bt4_model() + x0 = [ + 4.0382, + -2.9470, + -0.09115, + ] + + model = JuMP.Model() + @variable(model, x[i=1:3], start=x0[i]) + @objective(model, Min, x[1] - x[2] + x[2]^3) + @expression(model, circ, -25 + x[1]^2 + x[2]^2 + x[3]^2) + @constraint(model, circ == 0) + @constraint(model, x[1] + x[2] + x[3] == 1) + # Degenerate constraint + @constraint(model, circ - circ^2 == 0) + return MathOptNLPModel(model) +end + +function degen_30303_model() + lb = [0.0, 0.0, -Inf] + model = Model() + @variable(model, x[i=1:3] >= lb[i]) + @objective(model, Min, 0.1*x[1] + 0.1*x[2] + 0.8*x[3]) + @constraint(model, -x[1] - x[2] - x[3] + 1 >= 0) + @constraint(model, -x[1] - x[2] + x[3] + 1 >= 0) + @constraint(model, -x[1]*x[2] - (x[1] + x[2] + x[3] - 1)*(x[1] + x[2] - x[3] - 1) >= 0) + return MathOptNLPModel(model) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 1d3310c..004dd55 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,8 +18,21 @@ include("instances.jl") @test chk.inf_cc_con <= 1e-8 end +@testset "Test LICQ verification" begin + nlp = degenerate_bt4_model() + results = madnlp(nlp; print_level=MadNLP.ERROR) + degen_cons_bt4 = MadCheck.test_LICQ(nlp, results, MadCheck.ApproximatePrimalDualActiveSetLPEC(HiGHS.Optimizer), MadCheck.DegenJacSVD()) + + nlp = degen_30303_model() + results = madnlp(nlp; print_level=MadNLP.ERROR) + degen_cons_30303 = MadCheck.test_LICQ(nlp, results, MadCheck.ApproximatePrimalDualActiveSetLPEC(HiGHS.Optimizer), MadCheck.DegenJacSVD()) + + @test degen_cons_bt4 == [(constraints = [2, 3], bounds = [])] + @test degen_cons_30303 == [(constraints = [2, 3], bounds = [])] +end + @testset "Test active set detection" begin - nlp = hs15_model() + nlp = hs15_model() #TODO Find a better test to verify the methods are working properly - as they are inexact by nature, more over, the primal method can not detect weakly active constraints results = madnlp(nlp; print_level=MadNLP.ERROR) n = NLPModels.get_nvar(nlp) @@ -31,11 +44,11 @@ end new_point.solution[i] += (max_dist*noise)/n end - active_solution = MadCheck.find_active(nlp, results,MadCheck.ActiveMethodSimple()) + active_solution = MadCheck.find_active(nlp, results,MadCheck.BasicActiveSet()) - active_LP_P = MadCheck.find_active(nlp, new_point, MadCheck.ActiveMethodLP_P(HiGHS.Optimizer; Δ = 4 * 1e-8/n)) - active_LPEC = MadCheck.find_active(nlp, new_point, MadCheck.ActiveMethodLPEC(HiGHS.Optimizer)) - active_LPEC_A = MadCheck.find_active(nlp, new_point, MadCheck.ActiveMethodLPEC_A(HiGHS.Optimizer)) + active_LP_P = MadCheck.find_active(nlp, new_point, MadCheck.PrimalActiveSetLP(HiGHS.Optimizer; Δ = 4 * 1e-8/n)) + active_LPEC = MadCheck.find_active(nlp, new_point, MadCheck.PrimalDualActiveSetLPEC(HiGHS.Optimizer)) + active_LPEC_A = MadCheck.find_active(nlp, new_point, MadCheck.ApproximatePrimalDualActiveSetLPEC(HiGHS.Optimizer)) @test Set(active_LP_P.active) == Set(active_solution.active) && Set(active_LP_P.active_boundary) == Set(active_solution.active_boundary) @test Set(active_LPEC.active) == Set(active_solution.active) && Set(active_LPEC.active_boundary) == Set(active_solution.active_boundary)