[Nonlinear.ReverseAD] BoundsError from user-code #2523
Description
See https://discourse.julialang.org/t/max-of-a-vector-as-an-objective-for-jump/116911/8?u=odow
I haven't debugged yet
julia> using JuMP
julia> using Ipopt
julia> using Juniper
julia> function main()
β = 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)
end
main (generic function with 1 method)
julia> main()
nl_solver : MathOptInterface.OptimizerWithAttributes(Ipopt.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("print_level") => 0])
feasibility_pump : false
log_levels : [:Options, :Table, :Info]
#Variables: 4005
#IntBinVar: 1001
Obj Sense: Min
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:95
[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!
@ ~/.julia/packages/JuMP/7rBNn/src/optimizer_interface.jl:546 [inlined]
[22] main()
@ Main ./REPL[127]:50
[23] top-level scope
@ REPL[128]:1
Some type information was truncated. Use `show(err)` to see complete types.