Optimizing using objective with an integral in it

I’m new to JuMP and I’d like to minimize a function that has an integral in it. Here’s an example f(y) = \int_0^y y dy - y. Of course, this is just f(y)=(1/2)y^2-y which is minimized at y=1, but this is just a minimal example as a step towards something more complicated.

I don’t know how to pass this objective f into the JuMP stuff. I believe I need to register f as a “user-defined nonlinear function”, but I don’t know how to do that. According to the docs, “in addition to this list of functions, it is possible to register custom user-defined nonlinear functions.” But I don’t see how to do that and it doesn’t appear to say anything else.

Here’s the code I tried so far. The error message was: Unexpected object #89 (of type var"#89#90" in nonlinear expression.

Thanks!

using QuadGK
f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1] - y
using JuMP

function example()
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @NLobjective(model, Min, f)
    optimize!(model)

    return
end

example()
# ->  Unexpected object #89 (of type var"#89#90" in nonlinear expression.
1 Like

You should take a look at the JuMP documentation for some examples:

You need to define a variable @variable(model, y), you need to call the function f(y), and you’ll need to register the user-defined function register(model, :f, 1, f; autodiff = true).

I haven’t run the code so there might be a typo etc, but this should point you in the right direction:

function example()
    model = Model(Ipopt.Optimizer)
    register(model, :f, 1, f; autodiff = true)
    set_silent(model)
    @variable(model, y)
    @NLobjective(model, Min, f(y))
    optimize!(model)
    return
end

Depending on what you are trying to do, you might also want to take a look at GitHub - pulsipher/InfiniteOpt.jl: An intuitive modeling interface for infinite-dimensional optimization problems..

3 Likes

Thanks. From the docs here, I adjusted your code to replace the y's in your second line with foo and f. (using f instead of foo works too). This worked when I use the analytically defined function such as f(y) = 0.5*y^2 - y and gives a solution of y=1. But when I used the equivalent function (up to numerical approx) f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1] - y I get the following error saying something “repeats 79984 times”:

StackOverflowError:

Stacktrace:
 [1] cachedrule(#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1}}, n::Int64) (repeats 79984 times)
   @ QuadGK C:\Users\alex\.julia\packages\QuadGK\ENhXl\src\gausskronrod.jl:253

Here’s my code. Any thoughts on why it’s not working when I use the numerical integral version of the objective?

f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1] - y 
g(y) = 0.5*y^2 - y # g = f up to integration approx error
y = 0:0.1:10

function example()
    model = Model(Ipopt.Optimizer)
    @variable(model, y)
    register(model, :foo, 1, f; autodiff = true)
    set_silent(model)
    @NLobjective(model, Min, foo(y))
    optimize!(model)
    @show value(y)
    return
end

example()

Ha. That was an obvious typo. I’ve edited the post.

The issue is that your user-defined function needs to be differentiable:

julia> using ForwardDiff, QuadGK

julia> f(y) = quadgk(x -> x, 0, y, rtol=1e-8)[1]
f (generic function with 1 method)

julia> ForwardDiff.derivative(f, 0.0)
ERROR: StackOverflowError:
Stacktrace:
[1] cachedrule( #unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1}}, n::Int64 ) (repeats 79984 times)
@ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/gausskronrod.jl:253

Take a read of the posts here:

1 Like

QuadGK will not work like this with JuMP. You might try Quadrature.jl

1 Like

It appears that Quadrature.jl doesn’t yet support differentiation of the bounds of the integrals. Differentiation for a limit of an integral · Issue #56 · SciML/Quadrature.jl · GitHub

According to [Chirs R] Hasn’t been added yet. Someone would just need to extend Quadrature.jl/Quadrature.jl at master · SciML/Quadrature.jl · GitHub which wouldn’t be too hard.)