search

jump-dev / MathOptInterface.jl

A data structure for mathematical optimization problems
http://jump.dev/MathOptInterface.jl/
Other
387 stars 87 forks source link

Bug get/setting ConstraintSet with variable bridges #2452

Closed odow closed 5 months ago

odow commented 6 months ago

x-ref https://github.com/oxfordcontrol/Clarabel.jl/issues/160

julia> using JuMP, Clarabel

julia> model = Model(Clarabel.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Clarabel

julia> set_silent(model)

julia> @variable(model, x >= 1)
x

julia> @constraint(model, c, 2x == 3)
c : 2 x = 3

julia> optimize!(model)

julia> value(x)
1.4999999999999996

julia> b = normalized_rhs(c)
3.0

julia> set_normalized_rhs(c, b)

julia> optimize!(model)

julia> value(x)
0.9999999995620605
julia> import Clarabel

julia> import MathOptInterface as MOI

julia> src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
MOIU.UniversalFallback{MOIU.Model{Float64}}
fallback for MOIU.Model{Float64}

julia> x = MOI.add_variable(src)
MOI.VariableIndex(1)

julia> MOI.add_constraint(src, x, MOI.GreaterThan(1.0))
MathOptInterface.ConstraintIndex{MathOptInterface.VariableIndex, MathOptInterface.GreaterThan{Float64}}(1)

julia> c = MOI.add_constraint(src, 2.0 * x, MOI.EqualTo(3.0))
MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}(1)

julia> dest = MOI.instantiate(Clarabel.Optimizer; with_bridge_type = Float64)
MOIB.LazyBridgeOptimizer{MOIU.CachingOptimizer{Clarabel.MOIwrapper.Optimizer{Float64}, MOIU.UniversalFallback{MOIU.Model{Float64}}}}
with 0 variable bridges
with 0 constraint bridges
with 0 objective bridges
with inner model MOIU.CachingOptimizer{Clarabel.MOIwrapper.Optimizer{Float64}, MOIU.UniversalFallback{MOIU.Model{Float64}}}
  in state EMPTY_OPTIMIZER
  in mode AUTOMATIC
  with model cache MOIU.UniversalFallback{MOIU.Model{Float64}}
    fallback for MOIU.Model{Float64}
  with optimizer Empty Clarabel - Optimizer

julia> index_map = MOI.copy_to(dest, src)
MathOptInterface.Utilities.IndexMap with 3 entries:
  MOI.VariableIndex(1)                                           => MOI.VariableIndex(-1)
  ConstraintIndex{ScalarAffineFunction{Float64}, EqualTo{Float6… => ConstraintIndex{ScalarAffineFunction{Float64}, EqualTo{Float64}}(1)
  ConstraintIndex{VariableIndex, GreaterThan{Float64}}(1)        => ConstraintIndex{VariableIndex, GreaterThan{Float64}}(-1)

julia> MOI.set(dest, MOI.ConstraintSet(), index_map[c], MOI.EqualTo(3.0))

julia> MOI.optimize!(dest)
-------------------------------------------------------------
           Clarabel.jl v0.7.1  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 1
  constraints   = 2
  nnz(P)        = 0
  nnz(A)        = 2
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : Nonnegative = 1,  numel = 1

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, 
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0   0.0000e+00  -0.0000e+00  0.00e+00  4.46e-01  8.28e-01  1.00e+00  1.00e+00   ------   
  1   0.0000e+00   1.0946e+01  1.09e+01  2.85e-01  2.73e-02  1.14e+01  3.59e-02  9.66e-01  
  2   0.0000e+00   1.1412e+03  1.14e+03  2.85e-01  2.74e-04  1.14e+03  3.60e-04  9.90e-01  
  3   0.0000e+00   1.1417e+05  1.14e+05  2.85e-01  2.74e-06  1.14e+05  3.60e-06  9.90e-01  
  4   0.0000e+00   1.1417e+07  1.14e+07  2.85e-01  2.74e-08  1.14e+07  3.60e-08  9.90e-01  
  5   0.0000e+00   1.1417e+09  1.14e+09  2.85e-01  2.74e-10  1.14e+09  3.60e-10  9.90e-01  
---------------------------------------------------------------------------------------------
Terminated with status = primal infeasible
solve time = 2.41ms

julia> MOI.get(dest, MOI.VariablePrimal(), index_map[x])
0.9999999995620605
odow commented 6 months ago

Seems similar to https://github.com/jump-dev/MathOptInterface.jl/issues/2384 https://github.com/jump-dev/MathOptInterface.jl/pull/2396

odow commented 6 months ago

Here's a pure MOI

julia> import MathOptInterface as MOI

julia> src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
MOIU.UniversalFallback{MOIU.Model{Float64}}
fallback for MOIU.Model{Float64}

julia> x = MOI.add_variable(src)
MOI.VariableIndex(1)

julia> MOI.add_constraint(src, x, MOI.GreaterThan(1.0))
MathOptInterface.ConstraintIndex{MathOptInterface.VariableIndex, MathOptInterface.GreaterThan{Float64}}(1)

julia> c = MOI.add_constraint(src, 2.0 * x, MOI.EqualTo(3.0))
MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}(1)

julia> MOI.Utilities.@model(
           Model2452,
           (),
           (),
           (MOI.Nonnegatives, MOI.Zeros),
           (),
           (),
           (),
           (MOI.VectorOfVariables,),
           (MOI.VectorAffineFunction,)
       )
MathOptInterface.Utilities.GenericModel{T, MathOptInterface.Utilities.ObjectiveContainer{T}, MathOptInterface.Utilities.VariablesContainer{T}, Model2452FunctionConstraints{T}} where T

julia> function MOI.supports_constraint(
           ::Model2452{T},
           ::Type{MOI.VariableIndex},
           ::Type{
               <:Union{
                   MOI.GreaterThan{T},
                   MOI.LessThan{T},
                   MOI.EqualTo{T},
                   MOI.Interval{T},
                   MOI.ZeroOne,
                   MOI.Integer,
               },
           },
       ) where {T}
           return false
       end

julia> function MOI.supports_constraint(
           ::Model2452{T},
           ::Type{MOI.VectorOfVariables},
           ::Type{MOI.Reals},
       ) where {T}
           return false
       end

julia> dest = MOI.instantiate(Model2452{Float64}; with_bridge_type = Float64)
MOIB.LazyBridgeOptimizer{MOIU.GenericModel{Float64, MOIU.ObjectiveContainer{Float64}, MOIU.VariablesContainer{Float64}, Model2452FunctionConstraints{Float64}}}
with 0 variable bridges
with 0 constraint bridges
with 0 objective bridges
with inner model MOIU.GenericModel{Float64, MOIU.ObjectiveContainer{Float64}, MOIU.VariablesContainer{Float64}, Model2452FunctionConstraints{Float64}}

julia> index_map = MOI.copy_to(dest, src)
MathOptInterface.Utilities.IndexMap with 3 entries:
  MOI.VariableIndex(1)             => MOI.VariableIndex(-1)
  ConstraintIndex{VariableIndex, … => ConstraintIndex{VariableIndex, GreaterThan{Float64}}(-1)
  ConstraintIndex{ScalarAffineFun… => ConstraintIndex{ScalarAffineFunction{Float64}, EqualTo{Float64…

julia> MOI.get(dest, MOI.ConstraintSet(), index_map[c])
MathOptInterface.EqualTo{Float64}(1.0)

I can't help but think that many Variable bridges were a bad idea, particularly Vectorize ones with a non-zero constant.

odow commented 6 months ago

I think the issue is the interaction of the variable and constraint bridges, but it gets pretty complicated following the call_in_context.

mtanneau commented 6 months ago

Thanks for following up on that so quick!

What I was missing in trying to build a minimal example was that x should have a non-zero lower bound (I think, based on your example). Namely, I did not encounter any issue when x was declared as @variable(model, x >= 0). I think this echoes

particularly Vectorize ones with a non-zero constant

odow commented 6 months ago

A hacky work-around for now is remove_bridge:

julia> using JuMP, Clarabel

julia> model = Model(Clarabel.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Clarabel

julia> MOI.Bridges.remove_bridge(
           backend(model).optimizer,
           MOI.Bridges.Variable.VectorizeBridge{Float64},
       )

julia> set_silent(model)

julia> @variable(model, x >= 1)
x

julia> @constraint(model, c, 2x == 3)
c : 2 x = 3

julia> optimize!(model)

julia> value(x)
1.4999999999999996

julia> b = normalized_rhs(c)
3.0

julia> set_normalized_rhs(c, b)

julia> optimize!(model)

julia> value(x)
1.4999999999999996
odow commented 6 months ago

I need to discuss this with @blegat, but he's a bit busy this week.

odow commented 5 months ago

@blegat I think this one is for you. I really don't understand the ins and outs of call_in_context. I poke one bit and some unrelated stuff starts erroring.

odow commented 5 months ago

I swear I tested this:

julia> using JuMP, Clarabel

julia> model = Model(Clarabel.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Clarabel

julia> set_silent(model)

julia> @variable(model, x >= 1)
x

julia> @constraint(model, c, 2x == 3)
c : 2 x = 3

julia> optimize!(model)

julia> value(x)
1.4999999999999996

julia> b = normalized_rhs(c)
3.0

julia> set_normalized_rhs(c, b)

julia> optimize!(model)

julia> value(x)
0.9999999995620605

(clarabel) pkg> st
Status `/private/tmp/clarabel/Project.toml`
  [61c947e1] Clarabel v0.7.1
  [4076af6c] JuMP v1.21.0
  [b8f27783] MathOptInterface v1.27.1 `https://github.com/jump-dev/MathOptInterface.jl.git#master`
odow commented 5 months ago

Closed by #2472