From 86f96ac70249abc4dc30ca47bf582f6996bd83ed Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 21 Feb 2026 15:39:32 -0300 Subject: [PATCH] Cleanup Multiplicative Parameters handling --- src/MOI_wrapper.jl | 32 +++++++-------- src/ParametricOptInterface.jl | 3 -- src/cubic_objective.jl | 7 +++- src/duals.jl | 66 +++++++++++++++++++++++++++++++ test/test_JuMP.jl | 74 +++++++++++++++++++++++++++++++++++ test/test_MathOptInterface.jl | 12 +++++- test/test_cubic.jl | 73 ++++++++++++++++++++++++++++++++++ 7 files changed, 244 insertions(+), 23 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 0e605104..21e05b1f 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -74,10 +74,6 @@ function _cache_multiplicative_params!( for term in f.pv push!(model.multiplicative_parameters_pv, term.variable_1.value) end - for term in f.pp - push!(model.multiplicative_parameters_pp, term.variable_1.value) - push!(model.multiplicative_parameters_pp, term.variable_2.value) - end return end @@ -91,15 +87,22 @@ function _cache_multiplicative_params!( term.scalar_term.variable_1.value, ) end - for term in f.pp - push!( - model.multiplicative_parameters_pp, - term.scalar_term.variable_1.value, - ) - push!( - model.multiplicative_parameters_pp, - term.scalar_term.variable_2.value, - ) + return +end + +function _cache_multiplicative_params!( + model::Optimizer{T}, + f::ParametricCubicFunction{T}, +) where {T} + for term in f.pv + push!(model.multiplicative_parameters_pv, term.variable_1.value) + end + for term in f.pvv + push!(model.multiplicative_parameters_pv, term.index_1.value) + end + for term in f.ppv + push!(model.multiplicative_parameters_pv, term.index_1.value) + push!(model.multiplicative_parameters_pv, term.index_2.value) end return end @@ -137,7 +140,6 @@ function MOI.is_empty(model::Optimizer) isempty(model.vector_affine_constraint_cache) && # isempty(model.multiplicative_parameters_pv) && - isempty(model.multiplicative_parameters_pp) && isempty(model.dual_value_of_parameters) && model.number_of_parameters_in_model == 0 && isempty(model.parameters_in_conflict) && @@ -197,8 +199,6 @@ function MOI.empty!(model::Optimizer{T}) where {T} empty!(model.vector_affine_constraint_cache) # multiplicative_parameters_pv empty!(model.multiplicative_parameters_pv) - # multiplicative_parameters_pp - empty!(model.multiplicative_parameters_pp) # dual_value_of_parameters empty!(model.dual_value_of_parameters) # [SKIP] evaluate_duals diff --git a/src/ParametricOptInterface.jl b/src/ParametricOptInterface.jl index 5ade1c4c..08d2936d 100644 --- a/src/ParametricOptInterface.jl +++ b/src/ParametricOptInterface.jl @@ -186,7 +186,6 @@ mutable struct Optimizer{T,OT<:MOI.ModelLike} <: MOI.AbstractOptimizer ParametricVectorAffineFunction{T}, } multiplicative_parameters_pv::Set{Int64} - multiplicative_parameters_pp::Set{Int64} dual_value_of_parameters::Vector{T} evaluate_duals::Bool number_of_parameters_in_model::Int64 @@ -254,8 +253,6 @@ mutable struct Optimizer{T,OT<:MOI.ModelLike} <: MOI.AbstractOptimizer DoubleDicts.DoubleDict{ParametricVectorAffineFunction{T}}(), # multiplicative_parameters_pv Set{Int64}(), - # multiplicative_parameters_pp - Set{Int64}(), # dual_value_of_parameters T[], # evaluate_duals diff --git a/src/cubic_objective.jl b/src/cubic_objective.jl index 15bd0534..f9853777 100644 --- a/src/cubic_objective.jl +++ b/src/cubic_objective.jl @@ -45,10 +45,13 @@ function MOI.set( # 5. Clear old caches _empty_objective_function_caches!(model) - # 6. Store new cache + # 6. Cache multiplicative parameters + _cache_multiplicative_params!(model, cubic_func) + + # 7. Store new cache model.cubic_objective_cache = cubic_func - # 7. Store original for retrieval if option is enabled + # 8. Store original for retrieval if option is enabled if model.save_original_objective_and_constraints MOI.set( model.original_objective_cache, diff --git a/src/duals.jl b/src/duals.jl index 3bf8af47..eecc2d8b 100644 --- a/src/duals.jl +++ b/src/duals.jl @@ -16,6 +16,9 @@ function _compute_dual_of_parameters!(model::Optimizer{T}) where {T} if model.quadratic_objective_cache !== nothing _update_duals_from_objective!(model, model.quadratic_objective_cache) end + if model.cubic_objective_cache !== nothing + _update_duals_from_objective!(model, model.cubic_objective_cache) + end return end @@ -116,6 +119,69 @@ function _update_duals_from_objective!(model::Optimizer{T}, pf) where {T} return end +function _update_duals_from_objective!( + model::Optimizer{T}, + pf::ParametricQuadraticFunction{T}, +) where {T} + is_min = MOI.get(model.optimizer, MOI.ObjectiveSense()) == MOI.MIN_SENSE + sign = ifelse(is_min, one(T), -one(T)) + # p terms: ∂(c·p_i)/∂p_i = c + for term in pf.p + model.dual_value_of_parameters[p_val(term.variable)] += + sign * term.coefficient + end + # pp terms: ∂(c·p_i·p_j)/∂p_i = c·p_j + for term in pf.pp + mult = sign * term.coefficient + if term.variable_1 == term.variable_2 + mult /= 2 + end + model.dual_value_of_parameters[p_val(term.variable_1)] += + mult * model.parameters[p_idx(term.variable_2)] + model.dual_value_of_parameters[p_val(term.variable_2)] += + mult * model.parameters[p_idx(term.variable_1)] + end + return +end + +function _update_duals_from_objective!( + model::Optimizer{T}, + pf::ParametricCubicFunction{T}, +) where {T} + is_min = MOI.get(model.optimizer, MOI.ObjectiveSense()) == MOI.MIN_SENSE + sign = ifelse(is_min, one(T), -one(T)) + # p terms: ∂(c·p_i)/∂p_i = c + for term in pf.p + model.dual_value_of_parameters[p_val(term.variable)] += + sign * term.coefficient + end + # pp terms: ∂(c·p_i·p_j)/∂p_i = c·p_j (diagonal: c/2·2·p_i = c·p_i) + for term in pf.pp + mult = sign * term.coefficient + if term.variable_1 == term.variable_2 + mult /= 2 + end + model.dual_value_of_parameters[p_val(term.variable_1)] += + mult * model.parameters[p_idx(term.variable_2)] + model.dual_value_of_parameters[p_val(term.variable_2)] += + mult * model.parameters[p_idx(term.variable_1)] + end + # ppp terms: ∂(c·p_i·p_j·p_k)/∂p_i = c·p_j·p_k + for term in pf.ppp + coef = sign * term.coefficient + p1_val = model.parameters[p_idx(term.index_1)] + p2_val = model.parameters[p_idx(term.index_2)] + p3_val = model.parameters[p_idx(term.index_3)] + model.dual_value_of_parameters[p_val(term.index_1)] += + coef * p2_val * p3_val + model.dual_value_of_parameters[p_val(term.index_2)] += + coef * p1_val * p3_val + model.dual_value_of_parameters[p_val(term.index_3)] += + coef * p1_val * p2_val + end + return +end + function MOI.get( model::Optimizer{T}, attr::MOI.ConstraintDual, diff --git a/test/test_JuMP.jl b/test/test_JuMP.jl index abfb182c..b373596f 100644 --- a/test/test_JuMP.jl +++ b/test/test_JuMP.jl @@ -576,6 +576,80 @@ function test_jump_dual_objective_max() return end +function test_jump_dual_objective_pv() + # p*x is multiplicative → dual query should error + model = Model(() -> POI.Optimizer(HiGHS.Optimizer)) + set_silent(model) + @variable(model, x >= 0) + @variable(model, p in Parameter(-2.0)) + @objective(model, Min, p * x + x^2) + optimize!(model) + @test_throws ErrorException dual(ParameterRef(p)) + return +end + +function test_jump_dual_objective_pp_same() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer)) + set_silent(model) + @variable(model, x >= 1) + @variable(model, p in Parameter(3.0)) + @objective(model, Min, x + p^2) + optimize!(model) + # ∂(p²)/∂p = 2p = 6 + @test dual(ParameterRef(p)) ≈ 6.0 atol = 1e-4 + return +end + +function test_jump_dual_objective_pp_different() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer)) + set_silent(model) + @variable(model, x >= 1) + @variable(model, p1 in Parameter(2.0)) + @variable(model, p2 in Parameter(3.0)) + @objective(model, Min, x + p1 * p2) + optimize!(model) + # ∂(p1*p2)/∂p1 = p2 = 3, ∂(p1*p2)/∂p2 = p1 = 2 + @test dual(ParameterRef(p1)) ≈ 3.0 atol = 1e-4 + @test dual(ParameterRef(p2)) ≈ 2.0 atol = 1e-4 + return +end + +function test_jump_dual_objective_mixed_terms() + # min x + 2p + p^2 s.t. x >= p, p = 3 + # Optimal: x* = 3, obj = 3 + 6 + 9 = 18 + # ∂f/∂p from obj = 2 + 2p = 8 + # Constraint x - p >= 0: dual λ = 1, ∂g/∂p = -1, contribution = 1 + # Total = 8 + 1 = 9 + model = Model(() -> POI.Optimizer(HiGHS.Optimizer)) + set_silent(model) + @variable(model, x) + @variable(model, p in Parameter(3.0)) + @constraint(model, x >= p) + @objective(model, Min, x + 2 * p + p^2) + optimize!(model) + @test objective_value(model) ≈ 18.0 atol = 1e-4 + @test dual(ParameterRef(p)) ≈ 9.0 atol = 1e-4 + return +end + +function test_jump_dual_objective_parameter_update() + # min x + p^2 s.t. x >= 1, p = 2 then p = 3 + # ∂(p^2)/∂p = 2p + # p=2: dual = 4, p=3: dual = 6 + model = Model(() -> POI.Optimizer(HiGHS.Optimizer)) + set_silent(model) + @variable(model, x) + @variable(model, p in Parameter(2.0)) + @constraint(model, x >= 1) + @objective(model, Min, x + p^2) + optimize!(model) + @test dual(ParameterRef(p)) ≈ 4.0 atol = 1e-4 + set_parameter_value(p, 3.0) + optimize!(model) + @test dual(ParameterRef(p)) ≈ 6.0 atol = 1e-4 + return +end + function test_jump_dual_multiple_parameters_1() model = Model(() -> POI.Optimizer(HiGHS.Optimizer)) set_silent(model) diff --git a/test/test_MathOptInterface.jl b/test/test_MathOptInterface.jl index 61c6eeb1..cc9d0058 100644 --- a/test/test_MathOptInterface.jl +++ b/test/test_MathOptInterface.jl @@ -1451,8 +1451,16 @@ function test_qp_objective_parameter_times_parameter() 0.0, atol = ATOL, ) - @test MOI.get(optimizer, MOI.ConstraintDual(), cy) == 0.0 - @test MOI.get(optimizer, MOI.ConstraintDual(), cz) == 0.0 + @test isapprox( + MOI.get(optimizer, MOI.ConstraintDual(), cy), + 1.0, + atol = ATOL, + ) + @test isapprox( + MOI.get(optimizer, MOI.ConstraintDual(), cz), + 1.0, + atol = ATOL, + ) MOI.set(optimizer, MOI.ConstraintSet(), cy, MOI.Parameter(2.0)) MOI.optimize!(optimizer) @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 2.0, atol = ATOL) diff --git a/test/test_cubic.jl b/test/test_cubic.jl index 4da52914..e7aa5357 100644 --- a/test/test_cubic.jl +++ b/test/test_cubic.jl @@ -1623,6 +1623,79 @@ function test_jump_cubic_direct_model_ppp() return end +# ============================================================================ +# Contribution of objective parameters in duals +# ============================================================================ + +function test_cubic_dual_ppp_terms() + # min x + p^3 s.t. x >= 1, p = 2 + # Optimal: x = 1, obj = 1 + 8 = 9 + # ∂(p^3)/∂p = 3p^2 = 12 + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + @variable(model, x) + @variable(model, p in MOI.Parameter(2.0)) + @constraint(model, x >= 1) + @objective(model, Min, x + p^3) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test dual(ParameterRef(p)) ≈ 12.0 atol = ATOL + return +end + +function test_cubic_dual_ppv_terms() + # min p1*p2*x + x^2 s.t. x >= 0, p1 = 1, p2 = -2 + # ppv is multiplicative → dual query should error + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + @variable(model, x >= 0) + @variable(model, p1 in MOI.Parameter(1.0)) + @variable(model, p2 in MOI.Parameter(-2.0)) + @objective(model, Min, p1 * p2 * x + x^2) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 1.0 atol = ATOL + @test_throws ErrorException dual(ParameterRef(p1)) + @test_throws ErrorException dual(ParameterRef(p2)) + return +end + +function test_cubic_dual_pvv_terms() + # min p*x^2 - 3x s.t. 0 <= x <= 10, p = 1 + # pvv is multiplicative → dual query should error + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(1.0)) + @objective(model, Min, p * x^2 - 3 * x) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 1.5 atol = ATOL + @test_throws ErrorException dual(ParameterRef(p)) + return +end + +function test_cubic_dual_ppp_three_distinct() + # min x + p1*p2*p3 s.t. x >= 1, p1=2, p2=3, p3=4 + # ∂/∂p1 = p2*p3 = 12 + # ∂/∂p2 = p1*p3 = 8 + # ∂/∂p3 = p1*p2 = 6 + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + @variable(model, x) + @variable(model, p1 in MOI.Parameter(2.0)) + @variable(model, p2 in MOI.Parameter(3.0)) + @variable(model, p3 in MOI.Parameter(4.0)) + @constraint(model, x >= 1) + @objective(model, Min, x + p1 * p2 * p3) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test dual(ParameterRef(p1)) ≈ 12.0 atol = ATOL + @test dual(ParameterRef(p2)) ≈ 8.0 atol = ATOL + @test dual(ParameterRef(p3)) ≈ 6.0 atol = ATOL + return +end + end # module TestCubic.runtests()