Indicator constraints in IQP?

Specific Domains Optimization (Mathematical)
1

The problem I’m trying to solve has grown to the point where it requires quadratic constraints and/or objective. When the problem was linear, I had used indicator constraints successfully. After switching to Juniper or Alpine to solve the quadratic problem, I get an exception (“LoadError: Unable to use IndicatorToMILPBridge because element 2 in the function has a non-finite domain…”) during optimize! if I include an indicator constraint. Is this possibe? Am I missing something?

Minimal example:

using JuMP
using HiGHS
using Alpine
using Ipopt
using LinearAlgebra: dot

ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
highs = optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false)

model = Model(optimizer_with_attributes(
			Alpine.Optimizer,
			"nlp_solver" => ipopt,
			"mip_solver" => highs))

@variable(model, x[1:10,1:2], Bin)

values = [3.0, 4.0]

@objective(model, Max, dot(values[1].*x[:,1], values[2].*x[:,2]) )

@constraint(model, x[5,1] --> { sum(x[6:10,1]) <= 1 })

optimize!(model)

Exception:

ERROR: LoadError: Unable to use IndicatorToMILPBridge because element 2 in the function has a non-finite domain: 0.0 + 1.0 MOI.VariableIndex(6) + 1.0 MOI.VariableIndex(7) + 1.0 MOI.VariableIndex(8) + 1.0 MOI.VariableIndex(9) + 1.0 MOI.VariableIndex(10)
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] final_touch(bridge::MathOptInterface.Bridges.Constraint.IndicatorToMILPBridge{Float64, MathOptInterface.VectorAffineFunction{Float64}, MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.LessThan{Float64}}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{Alpine.Optimizer})
    @ MathOptInterface.Bridges.Constraint C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Bridges\Constraint\bridges\indicator_to_milp.jl:200
  [3] _final_touch(bridges::OrderedCollections.OrderedSet{MathOptInterface.Bridges.Constraint.IndicatorToMILPBridge{Float64, MathOptInterface.VectorAffineFunction{Float64}, MathOptInterface.ACTIVATE_ON_ONE, MathOptInterface.LessThan{Float64}}}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{Alpine.Optimizer})
    @ MathOptInterface.Bridges.Constraint C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Bridges\Constraint\map.jl:326
  [4] final_touch(map::MathOptInterface.Bridges.Constraint.Map, model::MathOptInterface.Bridges.LazyBridgeOptimizer{Alpine.Optimizer})
    @ MathOptInterface.Bridges.Constraint C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Bridges\Constraint\map.jl:333
  [5] final_touch
    @ C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Bridges\bridge_optimizer.jl:461 [inlined]
  [6] final_touch
    @ C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Bridges\bridge_optimizer.jl:466 [inlined]
  [7] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{Alpine.Optimizer}, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}})
    @ MathOptInterface.Utilities C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Utilities\copy.jl:508
  [8] copy_to
    @ C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Bridges\bridge_optimizer.jl:453 [inlined]
  [9] optimize!
    @ C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\MathOptInterface.jl:84 [inlined]
 [10] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Alpine.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities C:\Users\jeff.emanuel\.julia\packages\MathOptInterface\IiXiU\src\Utilities\cachingoptimizer.jl:316
 [11] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JuMP C:\Users\jeff.emanuel\.julia\packages\JuMP\R53zo\src\optimizer_interface.jl:448
 [12] optimize!(model::Model)
    @ JuMP C:\Users\jeff.emanuel\.julia\packages\JuMP\R53zo\src\optimizer_interface.jl:409
 [13] top-level scope
    @ C:\dev\monarch-prodsched-service\cmd\prodsched\src\mwe.jl:23
in expression starting at C:\dev\monarch-prodsched-service\cmd\prodsched\src\mwe.jl:23

versioninfo() prints

Julia Version 1.9.2
Commit e4ee485e90 (2023-07-05 09:39 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 12 × Intel(R) Core(TM) i7-10850H CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 12 on 12 virtual cores
Environment:
  JULIA_NUM_THREADS = 12

Pkg.status:

  [fbe9abb3] AWS v1.90.3
  [07493b3f] Alpine v0.5.6
  [336ed68f] CSV v0.10.11
  [861a8166] Combinatorics v1.0.2
  [a93c6f00] DataFrames v1.6.1
⌃ [cd3eb016] HTTP v1.9.15
  [87dc4568] HiGHS v1.7.5
  [b6b21f68] Ipopt v1.5.1
  [682c06a0] JSON v0.21.4
  [0f8b85d8] JSON3 v1.13.2
  [4076af6c] JuMP v1.17.0
  [2ddba703] Juniper v0.9.1
  [b8f27783] MathOptInterface v1.23.0
⌃ [bac558e1] OrderedCollections v1.6.2
  [aea7be01] PrecompileTools v1.2.0
  [856f2bd8] StructTypes v1.10.0
⌃ [5c2747f8] URIs v1.5.0
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [37e2e46d] LinearAlgebra
  [8dfed614] Test
Info Packages marked with ⌃ have new versions available and may be upgradable.

Thanks for any guidance.

2

Hmm. Let me take a look. I perhaps would have expected this to work, but I can see why it might not because you’re trying to combine quite a few concepts, and Alpine might reformulate the problem in a way that isn’t compatible with our bridge reformulations.

Do you have access to Gurobi? It can solve IQP with indicator constraints.

Note that if you just have quadratic terms between binaries, they can be reformulated as linear:

model = Model()
@variable(model, x[1:10,1:2], Bin)
values = [3.0, 4.0]
# @objective(model, Max, dot(values[1].*x[:,1], values[2].*x[:,2]) )
@variable(model, y[1:10], Bin)
@constraint(model, y .<= x[:, 1])
@constraint(model, y .<= x[:, 2])
@constraint(model, y .>= x[:, 1] .+ x[:, 2] .- 1)
@objective(model, Max, sum(values[1] .* values[2] .* y))
1 Like
3

Thanks. I don’t have access to Gurobi. I’ll look at reformulating the quadratic components in my actual problem similarly.

1 Like
4

I have the same problem when adding indicators into my (nonlinear) problem. I’ve tried using Alpine as well as Juniper, and I also get the error:

ERROR: Unable to use IndicatorToMILPBridge because element 2 in the function has a non-finite domain: 0.0 + 1.0 MOI.VariableIndex(2003)
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] final_touch(bridge::MathOptInterface.Bridges.Constraint.IndicatorToMILPBridge{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…})
    @ MathOptInterface.Bridges.Constraint ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/Constraint/bridges/indicator_to_milp.jl:200
  [3] _final_touch(bridges::OrderedCollections.OrderedSet{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…})
    @ MathOptInterface.Bridges.Constraint ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/Constraint/map.jl:340
  [4] final_touch(map::MathOptInterface.Bridges.Constraint.Map, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…})
    @ MathOptInterface.Bridges.Constraint ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/Constraint/map.jl:347
  [5] final_touch
    @ ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/bridge_optimizer.jl:450 [inlined]
  [6] final_touch
    @ ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/bridge_optimizer.jl:455 [inlined]
  [7] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/aJZbq/src/Utilities/copy.jl:394
  [8] copy_to
    @ ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/bridge_optimizer.jl:442 [inlined]
  [9] optimize!
    @ ~/.julia/packages/MathOptInterface/aJZbq/src/MathOptInterface.jl:121 [inlined]
 [10] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/aJZbq/src/Utilities/cachingoptimizer.jl:321
 [11] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:595
 [12] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:546
 [13] top-level scope

My complete example is here:

using JuMP
using Ipopt
using Juniper
using Alpine
using HiGHS
using Plots

## Settings

β = 0.5 # infectivity rate
γ = 0.25 # recovery rate
υ_max = 0.5 # maximum intervention
v_max = 10.0 # maximum duration of an intervention
silent = false
t0 = 0.0 # start time
tf = 100.0 # final time
δt = 0.1 # timestep
T = Int(tf/δt) # number of timesteps
S₀ = 0.99 # initial susceptible population
I₀ = 0.01 # initial infected population

## Set up model

#ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level"=>0)
#optimizer = optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>ipopt)
#model = Model(optimizer)

ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
highs = optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false)
model = Model(
    optimizer_with_attributes(
        Alpine.Optimizer,
        "nlp_solver" => ipopt,
        "mip_solver" => highs,
    ),
)

## Declare variables

@variable(model, S[1:(T+1)])
@variable(model, I[1:(T+1)])
@variable(model, υ[1:(T+1)])
@variable(model, v[1:(T+1)], Bin)

## Initial conditions
@constraint(model, S[1]==S₀)
@constraint(model, I[1]==I₀)

## Constraints on variables
@constraint(model, [t=2:(T+1)], 0 ≤  S[t] ≤ 1)
@constraint(model, [t=2:(T+1)], 0 ≤  I[t] ≤ 1)

## Constraints on control parameters
@constraint(model, [t=1:(T+1)], v[t] --> {υ[t]==υ_max})
@constraint(model, [t=1:(T+1)], !v[t] --> {υ[t]==0})
@constraint(model, δt*sum(v) ≤ v_max)

## Define auxiliary variables for infection and recovery
@NLexpression(model, infection[t=1:T], (1-exp(-(1 - υ[t]) * β * I[t] * δt)) * S[t])
@NLexpression(model, recovery[t=1:T], (1-exp(-γ*δt)) * I[t]);

## Set up timesteps
@NLconstraint(model, [t=1:T], S[t+1] == S[t] - infection[t])
@NLconstraint(model, [t=1:T], I[t+1] == I[t] + infection[t] - recovery[t])

## Minimize maximum of infection by defining an auxiliary variable
@variable(model, I_peak)
@constraint(model, 0.0 ≤ I_peak ≤ 1.0)
@constraint(model, [t=1:(T+1)], I[t] ≤ I_peak)

## Objective function
@objective(model, Min, I_peak)

if silent
    set_silent(model)
end
optimize!(model)
5

The issue is that you have indicator constraints involving υ, which does not have bounds.

Try adding something reasonable like:

@variable(model, 0 <= υ[1:(T+1)] <= 100)
6

If I set bounds on my variables as you suggest, I get rid of the above error, but now have the following error:

ERROR: BoundsError: attempt to access 4005-element Vector{Float64} at index [1:6007]
Stacktrace:
  [1] throw_boundserror(A::Vector{Float64}, I::Tuple{UnitRange{Int64}})
    @ Base ./abstractarray.jl:737
  [2] checkbounds
    @ ./abstractarray.jl:702 [inlined]
  [3] _copyto_impl!(dest::Vector{Float64}, doffs::Int64, src::Vector{Float64}, soffs::Int64, n::Int64)
    @ Base ./array.jl:374
  [4] copyto!
    @ ./array.jl:368 [inlined]
  [5] copyto!
    @ ./array.jl:388 [inlined]
  [6] _reverse_mode(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, x::Vector{Float64})
    @ MathOptInterface.Nonlinear.ReverseAD ~/.julia/packages/MathOptInterface/aJZbq/src/Nonlinear/ReverseAD/reverse_mode.jl:57
  [7] eval_constraint_jacobian(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, J::SubArray{…}, x::Vector{…})
    @ MathOptInterface.Nonlinear.ReverseAD ~/.julia/packages/MathOptInterface/aJZbq/src/Nonlinear/ReverseAD/mathoptinterface_api.jl:226
  [8] eval_constraint_jacobian(evaluator::MathOptInterface.Nonlinear.Evaluator{…}, J::SubArray{…}, x::Vector{…})
    @ MathOptInterface.Nonlinear ~/.julia/packages/MathOptInterface/aJZbq/src/Nonlinear/evaluator.jl:165
  [9] eval_constraint_jacobian(model::Ipopt.Optimizer, values::Vector{Float64}, x::Vector{Float64})
    @ Ipopt ~/.julia/packages/Ipopt/bqp63/src/MOI_wrapper.jl:750
 [10] (::Ipopt.var"#eval_jac_g_cb#6"{…})(x::Vector{…}, rows::Vector{…}, cols::Vector{…}, values::Vector{…})
    @ Ipopt ~/.julia/packages/Ipopt/bqp63/src/MOI_wrapper.jl:828
 [11] _Eval_Jac_G_CB(n::Int32, x_ptr::Ptr{…}, ::Int32, ::Int32, nele_jac::Int32, iRow::Ptr{…}, jCol::Ptr{…}, values_ptr::Ptr{…}, user_data::Ptr{…})
    @ Ipopt ~/.julia/packages/Ipopt/bqp63/src/C_wrapper.jl:0
 [12] IpoptSolve(prob::IpoptProblem)
    @ Ipopt ~/.julia/packages/Ipopt/bqp63/src/C_wrapper.jl:442
 [13] optimize!(model::Ipopt.Optimizer)
    @ Ipopt ~/.julia/packages/Ipopt/bqp63/src/MOI_wrapper.jl:962
 [14] optimize!(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Ipopt.Optimizer})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/bridge_optimizer.jl:367
 [15] solve_root_incumbent_model(jp::Juniper.JuniperProblem)
    @ Juniper ~/.julia/packages/Juniper/HBPrQ/src/model.jl:58
 [16] optimize!(model::Juniper.Optimizer)
    @ Juniper ~/.julia/packages/Juniper/HBPrQ/src/MOI_wrapper/MOI_wrapper.jl:301
 [17] optimize!
    @ ~/.julia/packages/MathOptInterface/aJZbq/src/Bridges/bridge_optimizer.jl:367 [inlined]
 [18] optimize!
    @ ~/.julia/packages/MathOptInterface/aJZbq/src/MathOptInterface.jl:122 [inlined]
 [19] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/aJZbq/src/Utilities/cachingoptimizer.jl:321
 [20] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:595
 [21] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:546

using JuMP
using Ipopt
using Juniper

β = 0.5 # infectivity rate
γ = 0.25 # recovery rate
υ_max = 0.5 # maximum intervention
v_max = 10.0 # maximum duration of an intervention
silent = false
t0 = 0.0 # start time
tf = 100.0 # final time
δt = 0.1 # timestep
T = Int(tf/δt) # number of timesteps
S₀ = 0.99 # initial susceptible population
I₀ = 0.01 # initial infected population

## Set up model

ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level"=>0)
optimizer = optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>ipopt)
model = Model(optimizer)


## Declare variables

@variable(model, 0 ≤ S[1:(T+1)] ≤ 1)
@variable(model, 0 ≤ I[1:(T+1)] ≤ 1)
@variable(model, 0 ≤ υ[1:(T+1)] ≤ υ_max)
@variable(model, v[1:(T+1)], Bin)

## Initial conditions
@constraint(model, S[1]==S₀)
@constraint(model, I[1]==I₀)

## Constraints on control parameters
@constraint(model, [t=1:(T+1)], v[t] --> {υ[t] ≤  υ_max})
@constraint(model, [t=1:(T+1)], !v[t] --> {υ[t]==0})
@constraint(model, δt*sum(v) ≤ v_max)

## Define auxiliary variables for infection and recovery
@NLexpression(model, infection[t=1:T], (1-exp(-(1 - υ[t]) * β * I[t] * δt)) * S[t])
@NLexpression(model, recovery[t=1:T], (1-exp(-γ*δt)) * I[t]);

## Set up timesteps
@NLconstraint(model, [t=1:T], S[t+1] == S[t] - infection[t])
@NLconstraint(model, [t=1:T], I[t+1] == I[t] + infection[t] - recovery[t])

## Minimize maximum of infection by defining an auxiliary variable
@variable(model, 0 ≤ I_peak ≤ 1)
@constraint(model, [t=1:(T+1)], I[t] ≤ I_peak)

## Objective function
@objective(model, Min, I_peak)
optimize!(model)
7

Let’s keep the discussion of this model to a single thread: Max of a vector as an objective for JuMP? - #9 by odow