From c7321e6815a98b3510fb66f73362dd8eeddf95ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 12 Oct 2023 22:23:52 +0200 Subject: [PATCH 1/5] Refactor KKT --- src/KKT/KKT.jl | 69 +--------------------------------------------- src/model.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+), 68 deletions(-) diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl index c973b24..adb0119 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -67,75 +67,12 @@ function PolyJuMP._invalidate!(model::Optimizer) return end -function _add_to_system(system, lagrangian, ::SS.FullSpace, ::Bool) - return lagrangian -end - -function _add_to_system( - system, - lagrangian, - set::SS.AlgebraicSet, - maximization::Bool, -) - n = SS.nequalities(set) - if iszero(n) - return - end - DynamicPolynomials.@polyvar λ[1:n] - for i in eachindex(λ) - p = SS.equalities(set)[i] - SS.add_equality!(system, p) - if maximization - lagrangian = MA.add_mul!!(lagrangian, λ[i], p) - else - lagrangian = MA.sub_mul!!(lagrangian, λ[i], p) - end - end - return lagrangian -end - -function _add_to_system( - system, - lagrangian, - set::SS.BasicSemialgebraicSet, - maximization::Bool, -) - lagrangian = _add_to_system(system, lagrangian, set.V, maximization) - DynamicPolynomials.@polyvar σ[1:PolyJuMP._nineq(set)] - for i in eachindex(σ) - p = SS.inequalities(set)[i] - SS.add_equality!(system, σ[i] * p) - if maximization - lagrangian = MA.add_mul!!(lagrangian, σ[i]^2, p) - else - lagrangian = MA.sub_mul!!(lagrangian, σ[i]^2, p) - end - end - return lagrangian -end - function _square(x::Vector{T}, n) where {T} return T[(i + n in eachindex(x)) ? x[i] : x[i]^2 for i in eachindex(x)] end function PolyJuMP._optimize!(model::Optimizer{T}) where {T} - if isnothing(model.options.solver) - system = SS.AlgebraicSet{T,PolyJuMP.PolyType{T}}() - else - I = SS.PolynomialIdeal{T,PolyJuMP.PolyType{T}}() - system = SS.AlgebraicSet(I, model.options.solver) - end - if model.model.objective_sense == MOI.FEASIBILITY_SENSE - lagrangian = MA.Zero() - else - lagrangian = MA.mutable_copy(model.model.objective_function) - end - lagrangian = _add_to_system( - system, - lagrangian, - model.model.set, - model.model.objective_sense == MOI.MAX_SENSE, - ) + lagrangian, system = PolyJuMP.lagrangian_kkt(model.model, model.options.solver) x = MP.variables(model.model) if lagrangian isa MA.Zero model.solutions = [ @@ -149,10 +86,6 @@ function PolyJuMP._optimize!(model::Optimizer{T}) where {T} model.raw_status = "Lagrangian function is zero so any solution is optimal even if the solver reports a unique solution `0`." return end - ∇x = MP.differentiate(lagrangian, x) - for p in ∇x - SS.add_equality!(system, p) - end solutions = nothing try # We could check `SS.is_zero_dimensional(system)` but that would only work for Gröbner basis based solutions = collect(system) diff --git a/src/model.jl b/src/model.jl index 8d9ec70..5cb14ca 100644 --- a/src/model.jl +++ b/src/model.jl @@ -335,6 +335,7 @@ function MOI.supports( ) where {T} return true end + function MOI.set( model::Model{T}, ::MOI.ObjectiveFunction, @@ -343,3 +344,77 @@ function MOI.set( model.objective_function = _polynomial(model.variables, func) return end + +function _add_to_system(_, lagrangian, ::SS.FullSpace, ::Bool) + return lagrangian +end + +function _add_to_system( + system, + lagrangian, + set::SS.AlgebraicSet, + maximization::Bool, +) + n = SS.nequalities(set) + if iszero(n) + return + end + DynamicPolynomials.@polyvar λ[1:n] + for i in eachindex(λ) + p = SS.equalities(set)[i] + SS.add_equality!(system, p) + if maximization + lagrangian = MA.add_mul!!(lagrangian, λ[i], p) + else + lagrangian = MA.sub_mul!!(lagrangian, λ[i], p) + end + end + return lagrangian +end + +function _add_to_system( + system, + lagrangian, + set::SS.BasicSemialgebraicSet, + maximization::Bool, +) + lagrangian = _add_to_system(system, lagrangian, set.V, maximization) + DynamicPolynomials.@polyvar σ[1:PolyJuMP._nineq(set)] + for i in eachindex(σ) + p = SS.inequalities(set)[i] + SS.add_equality!(system, σ[i] * p) + if maximization + lagrangian = MA.add_mul!!(lagrangian, σ[i]^2, p) + else + lagrangian = MA.sub_mul!!(lagrangian, σ[i]^2, p) + end + end + return lagrangian +end + +function lagrangian_kkt(model::Model{T}, solver = nothing) where {T} + if isnothing(solver) + system = SS.AlgebraicSet{T,PolyJuMP.PolyType{T}}() + else + I = SS.PolynomialIdeal{T,PolyJuMP.PolyType{T}}() + system = SS.AlgebraicSet(I, solver) + end + if model.objective_sense == MOI.FEASIBILITY_SENSE + lagrangian = MA.Zero() + else + lagrangian = MA.mutable_copy(model.objective_function) + end + lagrangian = _add_to_system( + system, + lagrangian, + model.set, + model.objective_sense == MOI.MAX_SENSE, + ) + if !(lagrangian isa MA.Zero) + ∇x = MP.differentiate(lagrangian, MP.variables(model)) + for p in ∇x + SS.add_equality!(system, p) + end + end + return lagrangian, system +end From 6595df6504b33e5234a831586058880d0c3009e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 12 Oct 2023 22:32:24 +0200 Subject: [PATCH 2/5] Fix format --- src/KKT/KKT.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl index adb0119..7c07b6a 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -72,7 +72,8 @@ function _square(x::Vector{T}, n) where {T} end function PolyJuMP._optimize!(model::Optimizer{T}) where {T} - lagrangian, system = PolyJuMP.lagrangian_kkt(model.model, model.options.solver) + lagrangian, system = + PolyJuMP.lagrangian_kkt(model.model, model.options.solver) x = MP.variables(model.model) if lagrangian isa MA.Zero model.solutions = [ From 9ce303d4c7929156d6964f6046af655e6c391ce2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 20 Oct 2023 11:22:02 +0200 Subject: [PATCH 3/5] Refactor --- src/model.jl | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/model.jl b/src/model.jl index 5cb14ca..4f04864 100644 --- a/src/model.jl +++ b/src/model.jl @@ -392,29 +392,45 @@ function _add_to_system( return lagrangian end -function lagrangian_kkt(model::Model{T}, solver = nothing) where {T} +function lagrangian_kkt( + objective_sense::MOI.ObjectiveSense, + objective_function::MP.AbstractPolynomialLike{T}, + set; + solver = nothing, + variables = nothing, +) where {T} if isnothing(solver) system = SS.AlgebraicSet{T,PolyJuMP.PolyType{T}}() else I = SS.PolynomialIdeal{T,PolyJuMP.PolyType{T}}() system = SS.AlgebraicSet(I, solver) end - if model.objective_sense == MOI.FEASIBILITY_SENSE + if objective_sense == MOI.FEASIBILITY_SENSE lagrangian = MA.Zero() else - lagrangian = MA.mutable_copy(model.objective_function) + lagrangian = MA.mutable_copy(objective_function) end lagrangian = _add_to_system( system, lagrangian, - model.set, - model.objective_sense == MOI.MAX_SENSE, + set, + objective_sense == MOI.MAX_SENSE, ) if !(lagrangian isa MA.Zero) - ∇x = MP.differentiate(lagrangian, MP.variables(model)) + ∇x = MP.differentiate(lagrangian, variables) for p in ∇x SS.add_equality!(system, p) end end return lagrangian, system end + +function lagrangian_kkt(model::Model{T}, solver = nothing) where {T} + return lagrangian_kkt( + model.objective_sense, + model.objective_function, + model.set; + solver, + variables = MP.variables(model), + ) +end From ba1aa6764704c3661ce331cc2e1cd9ab3770ad5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 20 Oct 2023 17:11:38 +0200 Subject: [PATCH 4/5] Fixes --- src/model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/model.jl b/src/model.jl index 4f04864..61d884d 100644 --- a/src/model.jl +++ b/src/model.jl @@ -393,7 +393,7 @@ function _add_to_system( end function lagrangian_kkt( - objective_sense::MOI.ObjectiveSense, + objective_sense::MOI.OptimizationSense, objective_function::MP.AbstractPolynomialLike{T}, set; solver = nothing, @@ -417,7 +417,7 @@ function lagrangian_kkt( objective_sense == MOI.MAX_SENSE, ) if !(lagrangian isa MA.Zero) - ∇x = MP.differentiate(lagrangian, variables) + ∇x = MP.differentiate(lagrangian, MP.variables(lagrangian)) for p in ∇x SS.add_equality!(system, p) end From 09efb0d16752dfb78ca21e934936fd282ae081a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Nov 2026 16:46:27 +0100 Subject: [PATCH 5/5] Fix lagrangian_kkt without objective and use model variables The objective function of the model is `nothing` when no objective is set, which did not match `objective_function::MP.AbstractPolynomialLike`. The `variables` keyword was ignored and the gradient was taken with respect to the variables of the Lagrangian, which misses model variables not appearing in it. It is now a required keyword. --- src/model.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/model.jl b/src/model.jl index 61d884d..546fe55 100644 --- a/src/model.jl +++ b/src/model.jl @@ -397,7 +397,7 @@ function lagrangian_kkt( objective_function::MP.AbstractPolynomialLike{T}, set; solver = nothing, - variables = nothing, + variables, ) where {T} if isnothing(solver) system = SS.AlgebraicSet{T,PolyJuMP.PolyType{T}}() @@ -417,7 +417,7 @@ function lagrangian_kkt( objective_sense == MOI.MAX_SENSE, ) if !(lagrangian isa MA.Zero) - ∇x = MP.differentiate(lagrangian, MP.variables(lagrangian)) + ∇x = MP.differentiate(lagrangian, variables) for p in ∇x SS.add_equality!(system, p) end @@ -428,7 +428,7 @@ end function lagrangian_kkt(model::Model{T}, solver = nothing) where {T} return lagrangian_kkt( model.objective_sense, - model.objective_function, + something(model.objective_function, zero(PolyType{T})), model.set; solver, variables = MP.variables(model),