Max of a vector as an objective for JuMP?

Specific Domains Optimization (Mathematical)
9

I can reproduce, so I’ll take a look: [Nonlinear.ReverseAD] BoundsError from user-code · Issue #2523 · jump-dev/MathOptInterface.jl · GitHub

Here’s how I would write your model though:

using JuMP
using Ipopt
using Juniper

β = 0.5
γ = 0.25
υ_max = 0.5
v_max = 10.0
silent = false
t0 = 0.0
tf = 100.0
δt = 0.1
T = round(Int, tf / δt)
S₀ = 0.99
I₀ = 0.01

ipopt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
optimizer = optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt)
model = Model(optimizer)
@variables(model, begin
    0 <= S[1:(T+1)] <= 1
    0 <= I[1:(T+1)] <= 1
    0 <= υ[1:(T+1)] <= υ_max
    v[1:(T+1)], Bin
    0 <= I_peak <= 1
end)
fix(S[1], S₀)
fix(I[1], I₀)
@expressions(model, begin
    infection[t in 1:T], (1 - exp(-(1 - υ[t]) * β * I[t] * δt)) * S[t]
    recovery[t in 1:T], (1 - exp(-γ * δt)) * I[t]
end)
@constraints(model, begin
    [t in 1:T+1], υ[t] <= υ_max * v[t]
    δt * sum(v) <= v_max
    [t in 1:T], S[t+1] == S[t] - infection[t]
    [t in 1:T], I[t+1] == I[t] + infection[t] - recovery[t]
    [t in 1:T+1], I[t] <= I_peak
end)
@objective(model, Min, I_peak)
optimize!(model)
1 Like