jump-dev / MathOptInterface.jl

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

[Bridges] Method error with quadratic indicator constraints #2506

Closed LebedevRI closed 2 months ago

LebedevRI commented 2 months ago

I've trimmed the reproducer somewhat, but it should be reducible further.

using Combinatorics
using JuMP
import HiGHS
#import MiniZinc
import LinearAlgebra

JuMP._CONSTRAINT_LIMIT_FOR_PRINTING[] = 10

function main()
    NUM_SWITCHES = 3
    NUM_KEYFRAMES = 3

    NUM_FRAMES_BETWEEN_KEYFRAMES = 10 # ???
    NUM_INTERKEYFRAME_ARCS = NUM_KEYFRAMES - 1
    NUM_FRAMES = NUM_KEYFRAMES + NUM_FRAMES_BETWEEN_KEYFRAMES * NUM_INTERKEYFRAME_ARCS

    KEYFRAME_INDEXES = collect(1:1+NUM_FRAMES_BETWEEN_KEYFRAMES:NUM_FRAMES)
    @assert length(KEYFRAME_INDEXES) == NUM_KEYFRAMES

    model = Model(HiGHS.Optimizer);
    #model = Model(() -> MiniZinc.Optimizer{Float64}("highs"))

    set_time_limit_sec(model, 1.0)

    @variable(model, KeyFrameMask[1:NUM_FRAMES], Bin)
    @variable(model, TransitionalFrameMask[1:NUM_FRAMES], Bin)
    @variable(model, DummyFrameMask[1:NUM_FRAMES], Bin)
    @constraints(model, begin
        # Each frame is in one of three states
        [f=1:NUM_FRAMES], KeyFrameMask[f]+TransitionalFrameMask[f]+DummyFrameMask[f] == 1
        # We know the number of keyframes
        sum(KeyFrameMask) == NUM_KEYFRAMES
        # Let's hardcode which frames are keyframes
        [f=1:NUM_FRAMES], KeyFrameMask[f] == (f in KEYFRAME_INDEXES)
        # In each inter-frame ark between keyframes, dummy frames are grouped right before next keyframe
        [g=1:NUM_INTERKEYFRAME_ARCS, f=1+KEYFRAME_INDEXES[g]+1:KEYFRAME_INDEXES[g+1]-1], DummyFrameMask[f-1] <= DummyFrameMask[f]
    end);    

    @variable(model, SwitchDrivePower[1:NUM_FRAMES, 1:NUM_SWITCHES], Bin)
    @constraints(model, begin
        # At most one switch can be operated per frame
        [f=1:NUM_FRAMES], 0 <= sum(SwitchDrivePower[f,:]) <= 1
        # At keyframes, switches are not operated
        [f=1:NUM_FRAMES, s=1:NUM_SWITCHES], KeyFrameMask[f] --> { SwitchDrivePower[f,s] == 0 }
        # At dummyframes, switches are not operated
        [f=1:NUM_FRAMES, s=1:NUM_SWITCHES], DummyFrameMask[f] --> { SwitchDrivePower[f,s] == 0 }
        # At transitionalframes, switches *are* operated
        [f=1:NUM_FRAMES], sum(SwitchDrivePower[f,:]) == TransitionalFrameMask[f]
        ### IF this is instead:
            # [f=1:NUM_FRAMES], TransitionalFrameMask[f] --> { sum(SwitchDrivePower[f,:]) == 1 }
        ### then it works

    end);

    # -1 - neutral; 0 - reverse; 2 - forward
    @variable(model, -1 <= SwitchDriveState[1:NUM_FRAMES, 1:NUM_SWITCHES] <= 2, Int)
    @constraints(model, begin
        # If switch's drive is power-on, there *IS* some control input
        [f=1:NUM_FRAMES, s=1:NUM_SWITCHES],  SwitchDrivePower[f,s] --> { SwitchDriveState[f,s] >=  0 }
        # If switch's drive is power-off, there is also no control input
        [f=1:NUM_FRAMES, s=1:NUM_SWITCHES], !SwitchDrivePower[f,s] --> { SwitchDriveState[f,s] == -1 }
    end);

    @variable(model, -1 <= SwitchDriveDirection[1:NUM_FRAMES, 1:NUM_SWITCHES] <= 1, Int)
    @constraints(model, begin
        # If switch's drive is power-off, there is also no control input
        [f=1:NUM_FRAMES, s=1:NUM_SWITCHES], !SwitchDrivePower[f,s] --> { SwitchDriveDirection[f,s] == 0 }
        # If switch's drive is power-on, there *IS* some control input
        [f=1:NUM_FRAMES, s=1:NUM_SWITCHES],  SwitchDrivePower[f,s] --> { SwitchDriveDirection[f,s] == SwitchDrivePower[f,s] * SwitchDriveState[f,s] - 1 }
    end);

    print(model)
    optimize!(model)
    solution_summary(model)
    @assert is_solved_and_feasible(model)
    #@show convert.(Bool, value.(KeyFrameMask))
end

main()
MethodError: no method matching bridge_constraint(::Type{MathOptInterface.Bridges.Constraint.IndicatorSOS1Bridge{Float64, MathOptInterface.EqualTo{Float64}}}, ::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, ::MathOptInterface.VectorQuadraticFunction{Float64}, ::MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.EqualTo{Float64}})

Closest candidates are:
  bridge_constraint(::Type{MathOptInterface.Bridges.Constraint.SplitIntervalBridge{T, F, S, LS, US}}, ::MathOptInterface.ModelLike, ::F, ::S) where {T, F, S, LS, US}
   @ MathOptInterface ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/Constraint/bridges/interval.jl:94
  bridge_constraint(::Type{MathOptInterface.Bridges.Constraint.FunctionConversionBridge{T, F, G, S}}, ::MathOptInterface.ModelLike, ::G, ::S) where {T, F, G, S}
   @ MathOptInterface ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/Constraint/bridges/functionize.jl:209
  bridge_constraint(::Type{MathOptInterface.Bridges.Constraint.IndicatorSOS1Bridge{T, S}}, ::MathOptInterface.ModelLike, ::MathOptInterface.VectorAffineFunction{T}, ::MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, S}) where {T<:Real, S}
   @ MathOptInterface ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/Constraint/bridges/indicator_sos.jl:43
  ...

Stacktrace:
  [1] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, BridgeType::Type, f::MathOptInterface.VectorQuadraticFunction{Float64}, s::MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.EqualTo{Float64}})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:1786
  [2] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, f::MathOptInterface.VectorQuadraticFunction{Float64}, s::MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.EqualTo{Float64}})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:1916
  [3] _copy_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap, index_map_FS::MathOptInterface.Utilities.DoubleDicts.IndexDoubleDictInner{MathOptInterface.VectorQuadraticFunction{Float64}, MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.EqualTo{Float64}}}, cis_src::Vector{MathOptInterface.ConstraintIndex{MathOptInterface.VectorQuadraticFunction{Float64}, MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.EqualTo{Float64}}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:259
  [4] _copy_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap, cis_src::Vector{MathOptInterface.ConstraintIndex{MathOptInterface.VectorQuadraticFunction{Float64}, MathOptInterface.Indicator{MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.EqualTo{Float64}}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:271
  [5] pass_nonvariable_constraints_fallback(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap, constraint_types::Vector{Tuple{Type, Type}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:282
  [6] pass_nonvariable_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, idxmap::MathOptInterface.Utilities.IndexMap, constraint_types::Vector{Tuple{Type, Type}})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:445
  [7] _pass_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, index_map::MathOptInterface.Utilities.IndexMap, variable_constraints_not_added::Vector{Any})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:330
  [8] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:505
  [9] copy_to
    @ ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:455 [inlined]
 [10] optimize!
    @ ~/.julia/packages/MathOptInterface/2CULs/src/MathOptInterface.jl:84 [inlined]
 [11] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{HiGHS.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/cachingoptimizer.jl:316
 [12] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/Gwn88/src/optimizer_interface.jl:457
 [13] optimize!
    @ ~/.julia/packages/JuMP/Gwn88/src/optimizer_interface.jl:409 [inlined]
 [14] main()
    @ Main ./In[12]:74
 [15] top-level scope
    @ In[12]:80
odow commented 2 months ago

Reproducer:

julia> using JuMP, HiGHS

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

julia> @variable(model, x)
x

julia> @variable(model, z, Bin)
z

julia> @constraint(model, z --> {x * z == 0})
z --> {x*z = 0}

julia> optimize!(model)
Running HiGHS 1.7.0 (git hash: 50670fd4c): Copyright (c) 2024 HiGHS under MIT licence terms
ERROR: MethodError: no method matching bridge_constraint(::Type{…}, ::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, ::MathOptInterface.VectorQuadraticFunction{…}, ::MathOptInterface.Indicator{…})

Closest candidates are:
  bridge_constraint(::Type{MathOptInterface.Bridges.Constraint.VectorSlackBridge{T, F, S}}, ::Any, ::MathOptInterface.AbstractVectorFunction, ::S) where {T, F, S}
   @ MathOptInterface ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/Constraint/bridges/slack.jl:351
  bridge_constraint(::Type{<:MathOptInterface.Bridges.Constraint.MultiSetMapBridge{T, S1, G}}, ::MathOptInterface.ModelLike, ::G, ::S1) where {T, S1, G}
   @ MathOptInterface ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/Constraint/set_map.jl:37
  bridge_constraint(::Type{MathOptInterface.Bridges.Constraint.NumberConversionBridge{T, F1, S1, F2, S2}}, ::MathOptInterface.ModelLike, ::F1, ::S1) where {T, F1, S1, F2, S2}
   @ MathOptInterface ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/Constraint/bridges/number_conversion.jl:46
  ...

Stacktrace:
  [1] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, BridgeType::Type, f::MathOptInterface.VectorQuadraticFunction{…}, s::MathOptInterface.Indicator{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:1786
  [2] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorQuadraticFunction{…}, s::MathOptInterface.Indicator{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:1916
  [3] _copy_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, index_map_FS::MathOptInterface.Utilities.DoubleDicts.IndexDoubleDictInner{…}, cis_src::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:259
  [4] _copy_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, cis_src::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:271
  [5] pass_nonvariable_constraints_fallback(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, constraint_types::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:282
  [6] pass_nonvariable_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, idxmap::MathOptInterface.Utilities.IndexMap, constraint_types::Vector{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:445
  [7] _pass_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, variable_constraints_not_added::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:330
  [8] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/copy.jl:505
  [9] copy_to
    @ ~/.julia/packages/MathOptInterface/2CULs/src/Bridges/bridge_optimizer.jl:455 [inlined]
 [10] optimize!
    @ ~/.julia/packages/MathOptInterface/2CULs/src/MathOptInterface.jl:84 [inlined]
 [11] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/2CULs/src/Utilities/cachingoptimizer.jl:316
 [12] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/Gwn88/src/optimizer_interface.jl:457
 [13] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/Gwn88/src/optimizer_interface.jl:409
 [14] top-level scope
    @ REPL[264]:1
Some type information was truncated. Use `show(err)` to see complete types.
LebedevRI commented 2 months ago

That reproducer results in

Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
KERNEL EXCEPTION
MethodError: no method matching namemap(::Type{MathOptInterface.ActivationCondition})
The applicable method may be too new: running in world age 31466, while current world is 31501.

Closest candidates are:
  namemap(::Type{MathOptInterface.ActivationCondition}) (method too new to be called from this world context.)
   @ MathOptInterface Enums.jl:214
  namemap(::Type{MathOptInterface.OptimizationSense}) (method too new to be called from this world context.)
   @ MathOptInterface Enums.jl:214
  namemap(::Type{LibGit2.Consts.GIT_TRACE_LEVEL})
   @ LibGit2 Enums.jl:214
  ...

which seems like a different issue from the one i posted.

odow commented 2 months ago

No, it results in the same issue and it is the reproducer. I have identified the fix.

I don't know what you've done to your Julia instance, but you must be @evaling some code somewhere to get the world age issue.

odow commented 2 months ago

Fix is https://github.com/jump-dev/MathOptInterface.jl/pull/2507

Your underlying issue is that you cannot solve this constraint with HiGHS.jl. The next version of MOI will throw a better error instead of the uninformative MethodError.

[f=1:NUM_FRAMES, s=1:NUM_SWITCHES],  SwitchDrivePower[f,s] --> { SwitchDriveDirection[f,s] == SwitchDrivePower[f,s] * SwitchDriveState[f,s] - 1 }
LebedevRI commented 2 months ago

@odow thank you for taking a look!

Yup, i very much expected that it would be too good to work :) The more faithful repro, which also does not work is

using JuMP, HiGHS
model = Model(HiGHS.Optimizer)
@variable(model, -1 <= x <= 1, Int)
@variable(model, -1 <= y <= 1, Int)
@variable(model, z, Bin)

# Variant using indicator constraints, works, but indicator constraints...
#@constraint(model, !z --> { y == 0 } )
#@constraint(model, z --> { y == x - 1 } )

@constraint(model, y == z * (x-1))
# You'd think it'd work, but:
# Constraints of type MathOptInterface.ScalarQuadraticFunction{Float64}-in-MathOptInterface.EqualTo{Float64} are not supported by the solver.

optimize!(model)
odow commented 2 months ago

HiGHS does not support quadratic constraints.

LebedevRI commented 2 months ago

Yes, i've understood that from the docs. (I mean, one could argue that they should be supported via falling back to the non-linear constraints.) I'm just showing the true snippet i was going for. @odow thank you!

odow commented 2 months ago

I mean, one could argue that they should be supported via falling back to the non-linear constraints

HiGHS does not support nonlinear constraints either. It is a mixed-integer linear solver, with additional support for continuous linear problems with a quadratic objective function.

LebedevRI commented 2 months ago

@odow thank you.

Question: would it be useful to file an issue here about potentially bridging

using JuMP, HiGHS
model = Model(HiGHS.Optimizer)
@variable(model, -1 <= x <= 1, Int)
@variable(model, -1 <= y <= 1, Int)
@variable(model, z, Bin)

@constraint(model, y == z * x)
   =>
@constraint(model, !z --> { y == 0 } )
@constraint(model, z --> { y == x } )

?

I mean, one could argue that they should be supported via falling back to the non-linear constraints

HiGHS does not support nonlinear constraints either. It is a mixed-integer linear solver, with additional support for continuous linear problems with a quadratic objective function.

Again, i get that :) What i meant is, there seems to be (some?) support in JuMP for user-defined functions, so perhaps the unsupported function types could be wrapped into such functions instead.

odow commented 2 months ago

would it be useful to file an issue here about potentially bridging

No. We won't bridge constraints based on some numerical structure. Bridges must (with some limited exceptions) support any input data. In this case, you are asking if ScalarQuadraticFunction-in-EqualTo can be bridged to a collection of affine indicator constraints. This is only possible in your case where at least one element in each quadratic term is binary.

Instead of trying to write JuMP models using the high-level syntax like --> etc, if you know that you want to use HiGHS to solve the problem then I strongly encourage you to write out the MIP reformulation by hand.