From 82c38dfb2379e213ee41ac09266c0f8c262d9273 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/4] 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 6be16ab..30dd54e 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -125,75 +125,12 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) return MOI.Utilities.default_copy_to(dest, src) 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 _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 = [ @@ -207,10 +144,6 @@ function _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 96499c4..67ca0fe 100644 --- a/src/model.jl +++ b/src/model.jl @@ -208,6 +208,7 @@ function MOI.supports( ) where {T} return true end + function MOI.set( model::Model{T}, ::MOI.ObjectiveFunction, @@ -216,3 +217,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 820e7f5c25e1e1ff0606cab9feae659babc1500a 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/4] 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 30dd54e..b99e1a6 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -130,7 +130,8 @@ function _square(x::Vector{T}, n) where {T} end function _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 3d9a0a3bf6e1274c0ce2e53011fa064aaf206db2 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/4] Refactor --- src/model.jl | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/model.jl b/src/model.jl index 67ca0fe..d68041b 100644 --- a/src/model.jl +++ b/src/model.jl @@ -265,29 +265,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 8fa92b4f8d1f84d23448b0a45e7daabf2d4ee8ff 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/4] Fixes --- src/model.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/model.jl b/src/model.jl index d68041b..add5871 100644 --- a/src/model.jl +++ b/src/model.jl @@ -266,7 +266,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, @@ -290,7 +290,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